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ABSTRACT 


A linear control subroutine library was created and 
stored in a load module on Disk 02 of the IBM/360 of the 
Naval Postgraduate School. 

This library consists of three groups of programs: 
transfer function subprograms; matrix manipulation and time 
response subprograms; and modern control design routines. 
The transfer function subprograms provide numerical aids 
for classical control design techniques including root locus 
and frequency design methods. The matrix manipulation and 
time response routines allow the user to determine eigen- 
values, find state transition matrices, evaluate resolvent 
matrices, perform several other matrix operations and 
determine and plot graphical time responses. The modern 
control design programs aid in solving Linear Quadratic 
Gaussian (LQG) problems and also provide the capability to 
investigate sensitivity and to de-couple multi-input multi- 
output systems. 

This thesis is a user's manual for the library of control 
design programs. Applications, extensive documentation and 
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EES 


A large number of computer programs are available to 
help today's control system engineers analyse and design 
increasingly complex systems. Most of these programs, how- 
ever, are available in the form of listings only, so any 
one wishing to use them must punch the cards, compile and 
test the routines, modify them and, most of the time, load 
them everytime a problem is to be solved. Obviously this is 
not a very practical and efficient way of using computers. 

The intent of this thesis was to facilitate the use of 
several of these programs by making them easily accessible 
ESL users as a pre-compiled load module library. The 
features of the library were to be as follows: 

- easy access to the subprograms 

- only rudimentary knowledge of FORTRAN coding and card 

set up procedures required to use the subprograms 

- good documentation readily available to the 

users (complete with subprogram descriptions, card 
set up procedures and worked out examples) 

- an expandable and improvable library 

- good priority eategories (class A or B only) 

for quick turnaround. 

Using these features as guidelines, a linear control sub- 

routines library (LINCON) was created, tested and is now 


available to any user on the Naval Postgraduate School IBM/360. 
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The following chapters describe the computer procedures 
and present these linear control subprograms contained in 
the library. All the different aspects of deck preparation 
and job control cards are discussed. The most common error 
conditions that may occur while using the subprograms and 
the remedial actions to be taken are pointed out. The linear 
control subprograms are presented in a user-oriented fashion. 
They are first introduced by defining their purposes and 
indicating the general rules that apply. Then the subpro- 
grams are individually described. The input requirements 
and the output to be expected are presented in great detail. 
Several examples are worked out, complete with the control 
cards, the input data, the computer output and the interpre- 
beegions of the results. 

Note that the programming aspects of the work are not 
included in the presentation. Reference 1, the provenance 
of most of the subprograms that constitute the LINCON library, 
must be consulted in that regard, along with the actual 
listings of the subroutines. Also note that reference 1 
can be used as an alternative source of information in 
using the subprograns. 

However, Appendix A explains how the LINCON data sets 
were created and gives the job control cards required to 
modify, verify or erase the data sets. Information on how 
to recreate the library, should it become necessary, are given 
as well. Finally, Appendix B specifies the references from 


which the worked examples were taken. 
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II. COMPUTER PROCEDURE 


The programs and subprograms described in [1], or modified 
versions, together with a few locally written programs were 
assembled in a load module to form a subroutine library. A 
private user disk space was allocated on Disk02 of the Naval 
Postgraduate School IBM/360 computer to hold the partitioned 
data set and library procedures were defined and cataloged so 
the library could be accessed by any computer user under OS 
Batch. The details on how to access the library and use the 
subprograms and subroutines are presented in the following 
paragraphs and a complete description of the data set, along 
with pertinent computer procedure information, is given in 
Appendix A. 

The system was devised in such a way as to minimize the 
need for programming and provide the user with a convenient, 
flexible, easy-to-utilize tool for analysis and design of 
linear control systems. The programs and subprograms were 
kept as separate subroutines so one, or more, programs could 
be executed as a single job. The following gives a detailed 
description of the different methods of accessing the library 
as well as the cards necessary to run the programs under batch 
processina. For convenience, the major steps of the proce- 
dure are also reproduced in Section III as part of the sub- 


programs presentation. 
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A. MODES OF OPERATION AND CONTROL CARDS 

There are many different overating modes a user can 
employ to access a given computer library. Of these methods, 
three were determined to be most appropriate and are presented 
hereafter. It is pointed out that these procedures only 
apply to this specific set of programs, which was assigned the 
name LINCON (for linear control). Also observe that each mode 
implies the use of slightly different control card deck set 
ups. These differences are essential to the proper operation 
of the system under the selected mode of operation. Each line 
must be meticulously reproduced on the computer card and the 
order of appearance of the cards scrupulouslv respected. 

1. Mode One 

This mode applies when a user wants to execute only 
one of the subprograms for either single or multiple runs. 
Except for the subprograms named GTRESP, KALMAN and PRTLOC, 
which require exterior subroutines, all subprograms can be 
accessed using this method. Mode Two establishes the proce- 
dures that deal with the three special cases enumerated above. 
For Mode One, the control cards must be: 

// (standard OS JOB card) 
// EXEC, LINCON 
//LINK.SYSIN,DD * 
aa INCLUDE, SYSLIB (member) 
/* 


MWIGOMSYSTN, DD, * 
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data deck as described in Section III 

for the subprogram "member" 
pa^ 
where "member" is the simple name defining a subprogram to 
be executed. For example, "include syslib(SERCOM)" would 
have to be typed on the appropriate card to access the sub- 
routine library program called SERCOM. 

2. Mode Two 

The three special cases previously mentioned are 
accessed using this mode of operation. A different library 
procedure was created since GTRESP, KALMAN or PRTLOC might 
very well require special functions or inputs that vary as 
given parameters change. This situation does not significantly 
complicate the procedure and greatly adds to the system capa- 
bility. Further justification and explanation are given in 
Section III along with the subprogram descriptions. Again 
under this second mode, the programs are to be accessed one 
at a time, either for single or multiple runs. The computer 
card deck set up to be provided is: 
// (standard OS JOB card) 
// EXEC ,LINCONF 
//FORT.SYSIN,DD,* 

FORTRAN deck of user supplied subroutine as 

specified for GTRESP, KALMAN or PRTLOC 
/* 


EINE. SYSIN.DD _* 


Un 











an INCLUDE~SYSLIB (member) 

n ENTRY member 

/* 

77/GO.SYSIN.~DD.* 
data deck for "member" as described 
ENEE ep Lil. 

/* 

where member d the actual name of the subprogram 
to be executed. In this case, it is either GTRESP, KALMAN 
or PRTLOC. 

For example, "include syslib(KALMAN)" on the appro- 
priate card, followed by "entry KALMAN" on the next card 
would cause the subprogram called KALMAN to be run. 

3. Mode Three 

This mode of operation permits the user to call more 
than one subprogram while executing a single job. Since this 
third option calls all the subprograms at the same time, a 
large amount of computer memory is required. The user must 
be aware that this increases the turnaround time. Nevertheless 
the method can still be very useful. For instance, a user who 
is not in a hurry could utilize this set up to obtain the 
solution to several simple problems which do not require 
modification of some parameters. 

At this point the user is reminded that great care 
must be taken to correctly prepare the control and data decks. 
With an increased turnaround time, errors become costly and 


Very frustrating. 
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When it is decided to use Mode Three, the following 
computer cards must be generated: 
// (standard OS JOB card), TIME=5 
// XEXECALINCON,REGION.GO-2350K 


ES UNCLUDE .SYSLIB (MAIN) 


> 
//GO.SYSIN,DD,* 
MEMBER 1 
data deck for member 1 as 
described in Section III 
nr 2 
data deck for member 2 as 
described in Section III 
$ 
/* 


where MEMBER l, MEMBER 2, etc., are the defining names 
of the subprograms to be executed and start in column one. 
Note that again, as explained in Section III, the data deck 
pertaining to the same subprogram can be arranged either for 
single or multiple runs. The dollar sign, $, is a stop sign 
to be printed in column one. This dollar sign card must appear 
E tor the last data deck of each "member" to be executed 


under Mode Three. 


B. ERROR CONDITIONS 
When running programs, it is rather disappointing if 


results do not come out as expected. This in itself is a good 
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reason to always verify one last time that the control cards 
were punched correctly and the data deck was set up exactly 
as specified. Nonetheless, both neophytes and veterans do 
make mistakes and the purpose of this section is to outline 
some of the most common errors and show how to identify and 
correct them. The user must keep in mind that the error 
conditions and messages presented below apply to the IBM/360 
and were taken from [2] which is the only up-to-date source 
of information on the subject. 

Before any attempt is made to correct an eventual problem, 
the errors must be 'exposed'. This very important step is 
too often jumped over, the user opting to guess directly what 
went wrong. In order to save time and effort, one should pro- 
ceed more logically. The user should always check the linkage 
editor and job scheduler output to ascertain that the proper 
actions indeed did take place. Any messages such as '-Step- 
Go-Was Not Run Because of Condition Codes' clearly indicate 
what operations were not carried out and direct the user to 
the problem. Using these makes it much easier for the pro- 
grammer to pinpoint the malfunction and take the appropriate 
action. If no faulty indications appear in the messages out- 
put by the job scheduler (IEF type messages), the linkage 
editor (IEW type messages), the program producing (IEY) or 
the object program (IHC) and the results obtained are still 
suspected to be erroneous, the user then knows he should devote 


his attention to the mathematics of the problem and revise 
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Erezimpucr data (i.e., the output obtained is not the result 
era computer error'). 

Some of the possible linkage editor, object program and 
program Se messages are listed below. These should 
give the programmer a good idea of what to expect and how to 
proceed. Experience has shown that even if only a minimum of 
information is provided, the user greatly benefits from having 
these simple explanations at hand. 

l. IEWOOO (control statement only) 

This message enumerates all the control statements 
passed to the linkage editor. INCLUDE and ENTRY cards are 
listed for reference. it. is not an error message. 

2. IEWOl32 ERROR - SYMBOL PRINTED IN AN UNRESOLVED 

EXTERNAL REFERENCE 

This indicates that the symbol printed to the right 
of IEW0132 is a subprogram or subroutine which was not in the 
specified load module library or other modules passed to the 
linkage editor for processing. The user must make sure the 
correct subroutine library was specified (i.e., LINCON or 
LINCONF as required for proper mode of operation), and that 
the subroutine name requested was correctly spelled. 

3. 1EW0222 ERROR - CARD PRINTED CONTAINS INVALID INPUT 

FROM OBJECT MODULE. 

In this case, either some control cards were missing, 
thus causing the editor to interpret wrongly the cards that 
followed, or some of the cards were punched incorrectly. The 


deck should be checked. 
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4. IEW0342 - LIBRARY SPECIFIED DOES NOT CONTAIN MODULE. 
The subprogram or subroutine name specified on the 
INCLUDE control card was not found in the LINCON library. 
The user must make sure the INCLUDE card was punched as 
follows: 
INCLUDE SYSLIB (member) 
where 'member' is the appropriate subprogram name. 
5. IHC900I EXECUTION TERMINATING DUE TO ERROR COUNT FOR 
ERROR NUMBER 217 
MZA APETEOS — END OF DATA SET ON UNIT 5 
Here the computer stopped executing due to lack of 
data. At that instant, the problem might have been completely 
solved or not. It is advisable not to take any chances. Again 
the data deck should be thoroughly checked to ascertain that 
the cards were set up properly and the data deck incorporated 
was really the one for the specified subroutine. 
6. IHC215 CONVERT - ILLEGAL DECIMAL CHARACTER (decimal 
character) 
The computer found the given decimal character where 
a number was expected. Either the data cards were improperly 
set up, the subprogram name specified was incorrect or the 
FORTRAN format specified was not adhered to. Remedial actions 
Should be taken accordingly. 
7]. %<IEYO32I NUL PROGRAM 
This message indicates that no exterior subroutine 


was provided when needed and that the computer considered all 
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the data expected from this subroutine to be zero. Even if 
this situation can sometimes be used to advantage, it is not 
recommended here. The programmer should incorporate all 
required subroutines in his deck. Note that this error can 
only occur while accessing the library under Mode Two. 

8. No error condition messages printed out but incom- 
plete or no results were output by the computer. Here many 
things could have gone wrong, but most likely one of the 
bolllowing occurred: 

- While operating under Mode Two, the ENTRY card was not 
provided where required. The user must verify the program 
cards for correctness. 

- While operating under Mode Three, insufficient region 
size was specified. The remedial action is then to increase 
region size. 

- While operating under any of the three modes and the 
two conditions described above were not the cause, insuffi- 
cient running time was allocated for the program. If the 
CPU time indicated on the output and the one specified on 
the JOB card matched, the user should then allow more time 
for computation. 

- If none of the above, an underflow or overflow condi- 
tion may have occurred, causing the program to stop. In 
this case the linkage editor and job scheduler output will 
indicate a completion code - OCF. The user must verify the 


data cards and make the appropriate corrections. 
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The error conditions listed above are obviously not the 
only ones that can occur, but they are the ones a user is 
most likely to come across while employing the subroutine 


library called LINCON. 
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III. THE LINEAR CONTROL PROGRAMS 


pel INTRODUCTION 

The previous section dealt with the control statements 
and the card deck arrangements required to introduce the 
computer jobs to the operating system and tell the latter 
everything it needs to know about the input and output 
requirements. This chapter introduces the theory necessary to 
use the programs, presents a precise description of all sub- 
routines and data cards and gives detailed examples taken 
among problems that were solved on the computer. 

T. Outline 

The subprograms are divided into three classes: the 

transfer function subprograms, the time response and matrix 
manipulation subprograms and the modern control subprograms. 
The first set allows the user to obtain a root locus starting 
from a block diagram or signal flow graph (RTLOC), the 
Boots of a polynomial and their locus (PRTLOC), the 
Bode and Nyquist frequency plots (FRESP), the partial frac- 
tion expansion of the ratio of two polynomials (PRFEXP) and, 
finally, the roots of any polynomial (ROOTS). The second 
group is composed of three subprograms which are provided for 
determining the rational time response (RTRESP) and the 
graphical time response (GTRESP) of linear feedback control 
Systems and for computing the determinant, inverse, charac- 


teristic polynomial, eigenvalues, state transition matrix and 


25 





the resolvent matrix (BASMAT). The last group of subprograms 
deals with optimal control design. It permits the user to 
find the observability index of a control system (OBSERV), 
to test for both controllability and observability (CONOBS), 
to obtain the state variable feedback given some performance 
criterion (STVAR), to determine the complete sensitivity 
analysis of the closed-loop system poles variation as certain 
parameters change (SENSIT), to design Luenburger observers 
(LUEN) and serial compensators (SERCOM), to minimize a per- 
formance index when some state variables are inaccessible, 
to solve the Riccati equation to derive optimal control 
parameters and continuous Kalman filters (RICATI), to compute 
the gains of discrete Kalman filters (KALMAN), to evaluate 
the feedback control gains for discrete linear regulator 
problems, and, finally, to decouple multiple-input multiple- 
output systems (MIMO). Table I conveniently summarizes the 
above. 

The purpose of each subprogram and a brief discussion 
of the theory behind it are given in the subprograms presentation. 

2. Input Format 

The input format for each of the subprograms is com- 
pletely described with their presentation and must be referred 
to in each case. However, since the same general input for- 
mat is used for all the programs, it is appropriate to point 
out some of the similarities and the conventions adopted. 
For instance, to make it easier to remember, most of the 


groups of data cards have the same arrangement. 
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TABLE I - The Linear Control Subprograms 





Mode of 

Name Purpose ‘ Operation Class 

RTLOC To plot the root of a One or B 
polynomial equation Three 
starting from a feed- 
back control system 
block diagram. 

PRTLOC To plot the root locus Two | B 
of a characteristic d 
polynomial. 

FRESP To obtain and plot the One or A/B 
frequency response of a Three 
rational transfer 
function over a 
Specified range of 
frequencies. Both Bode 
and Nyquist diagrams can 
be plotted. 

PRFEXP TO perform the partial One or A 
fraction expansion of a Three 
rational function. 

ROOTS To find the roots of a One or A 


polynomial of order less Three 
than or equal to twenty. 





BASMAT To compute the determi- One or A 
nant, the inverse, the Three 
characteristic poly- 
nomial, the eigenvalues, 
the state transition 
matrix, and the 
resolvent matrix from 
a given matrix 
A (NxN). 





RTRESP To determine the rational One or A 
time response of a system Three 
(in closed-form). 
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TABLE I (Continued) 











Mode of 

Name Purpose Operation Class 

GTRESP To obtain the graphical One or A/B 
time response of a Three 
system for a specified 
inpūt. 

OBSERV To determine the One or A 
observability index for Three 
a linear system. 

CONOBS To check for both One or A 
observability and Three 
controllability of a 
linear system. 

SENSIT To study the closed- One or A/B 
loop poles variation Three 
of a linear feedback 
system. 

STVAR To calculate the One or A 
controller gain and Three 
feedback coefficients 
to achieve a desired 
closed-loop transfer 
function. Also computes 
the plant transfer 
function, internal 
transfer functions and 
determines EA (s), the 
equivalent single- 
feedback element. 

LUEN To design  Luenberger one or A 
Observers to achieve a Three 
desired closed-loop 
transfer function. 

SERCOM To deSign a series One or A 
compensator to achieve Three 


a desired closed-loop 
transfer function. 











TABLE I (Continued) 


Mode of 
Name Purpose Operation Class 
RICATI To solve the differential One or A/B 
Riccati equation to Three 
determine the optimum 
control gains for state- 
regulator problems and/or 
tne continuous Kalman 
filter gains. 
KALMAN To determine the discrete Two A 
Kalman filter gains. 
STREG To evaluate the discrete One or A 


feedback gains of linear Three 
regulator problems. 


E a a a a a m 


MIMO To decouple an Nth order One or A 
system with M inputs and Three 
M outputs and place the 
closed-loop poles of 
each decoupled sub- 
system at specified 
locations. 


nn nn 
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a. First Data Card 

The purpose of the first data card of any of the 
subprograms is to identify the problem for future reference 
and for output data. A maximum of twenty alpha-numeric 
characters (except $) can be used, starting in column one 
(format 5A4). On this first card, the user also normally 
defines the system order and the dimensions of the various 
matrices (format I2 for each number to be entered). Note 
that the dollar sign $ has been defined as a STOP and must 
never be used as problem identification. 

b. Matrices 

Matrices are entered one row at a time either in 
Menr original form or transposed, as specified. The input format 
table presented with each subprogram indicates the correct form to 
use. The vectors are always defined using lower case letters 
while other matrices are identified with capital letters. 
The matrix elements are punched in ten-column fields (format 
8El0 or 8F10), thus a maximum of eight numbers can be given 
per card. If the order of the system is greater than eight, 
two cards are needed for every row. 

An example will demonstrate the procedure. Assume 


that the A and b matrices are: 


^— be, d 


3.19 0.00 UE INS 
A = 25.45 6.40 = 0.50 5 - 0.0 
1.00 =9.14 6.75 1522 
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The given A and b matrices are entered using an 8F10.3 


format as follows: 


card columns a 11 ou 
3. 19 0.0 gomme 
A 2.45 6.4 ADO 
O -9,14 6:75 
b? 1.0 0.0 SC 


C. Polynomials 

The polynomial data can be entered in two differ- 
ent formats referred to as P mode (polynomial form) and F 
mode (factored form). If P mode is selected, the letter P 
(format Al) followed by the degree of the polynomial (format 
12) are entered on One card. The coefficients of the polynomial 
are placed on the next card(s) each in ten column fields 
(format 8F10 or 8E10). The polynomials are always presented 
in ascending order, the constant term given first and the 
coefficient of the highest term assumed to be unity. In other 
words, the last coefficient entered will always be interpreted 
as being 1.0, thus can be entered either as '1.0' or as a 
blank. Again an example best illustrates the principles. 

The given four polynomials are entered using an 
SmO = format: 

(1) Polynomials: 


TIO as s 


(11) st 53° 4 6s? + E 


DAL 





Ma). E, (highest degree coefficient 
of a zero order polynomial) 
(iv) 4 + T + at T 35° + e? 
(2) Computer data cards: 
card 1] Dun 21 31 41 51 61 pu 
columns 
P02 


If it is desired to enter the polynomial in 
factored form, then the F mode is chosen. This choice is 
indicated by placing the letter F (format Al) in the first 
column followed by the degree of the polynomial in the next 
two (format I2). The factors are then entered one per card, 
the real part in the first ten column field and the imaginary 
part in the next ten columns (format 2El0 or 2F10). An 
unusual convention was picked to enter all the possible factors. 
The user must be careful and make sure his notation agrees 
with the following: 

(1) The real part of the root is entered as 


positive if the factor is in the left half plane. 
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(2) The real part of the root is entered as 
negative if the factor is in the right half plane. 

(3) Only one of the complex conjugate roots is 
entered and it must be with the one with the positive 
imaginary part. 

(4) If the polynomial is a constant, it must 
equal 1.0 and be entered in the P mode. as shown below. 
Examples covering many possibilities are shown next. 

Factored polynomials: 

(i) (s + 3) 1) 

(ii) s(s + 4)(s + 1 + j)(s + 1 - j) 

Gain alao 

(iv) (s - 1l1)(s - 2 + j$)(s-2- j5) (s +3) (s +3) 

Computer data cards: | 


FO2 


221) 
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F04 


(iv) -2.0 520 


One good way to remember how to work around this 
confusing notation is to always enter the real parts as they 
appear in the factored polynomial and include the positive 
imaginary part only. In other words, one can analyse any 


situation in the following manner: 


EN Us 3) (s - 1) {s + 1+ 42.) (s + 1 - 32) 


where the circles indicate the numbers to be punched. 
d. Multiple Runs 

One last common characteristic of the input data 
is that one or several data decks pertaining to the same 
subprogram can be stacked and run as a single job. in other 
words, one complete data deck is prepared for each problem 
but the decks are all put one on top of the other and read in 
to the computer preceded only by one set of control cards. 

The user must realize,however,that this feature 
implies more runs to be performed in a single job and thus 
the time limit to be specified on the JOB control card must 


be estimated accordingly. 
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PEROU CPUE Format 


The output of all of the subprograms is quite com- 
prehensible and need not be explained. Nevertheless confusion 
may arise due to certain factors that are commented upon here. 
For the matrices, the same rules as for the input apply; 
vectors are listed out as transpose matrices and all other 
types of matrices are presented one row at a time. For 
convenience, polynomials are always output both in polynomial 
and factored forms no matter how they were provided as input. 
As for the input, the coefficients appear in ascending order, 
the constant term first. In factored form, the roots are 
listed with their normal sign convention; the left half 
plane roots are negative and those in the right half plane 
positive. 

Hence there is a sign inversion between the input 


and the output for the factored case. 


B. THE TRANSFER FUNCTION SUBPROGRAMS 
This set helps the user to analyse or design feedback 
control systems by providing a means of obtaining quickly the 
roots locus, Bode diagrams, Nyquist plots, partial fraction 
expansions and polynomial roots. 
1. Root Locus (RTLOC) 
This subprogram calculates and plots the roots of the 


equation 


l + K G(s) = 0 
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where G(s) is a rational function of the form 


Z 
D 





J 
u 


The user must provide N(s), D(s) and a range of value for 
K. Since a choice of two ways to vary K from minimum to 
maximum gain exists, an option card is also required. 
ac Input 
The observations and the table presented below 
should be sufficient to use the subproaram which can be 
called under Mode One or Mode Three (as described in Chapter 
E : 
(1) N(S) can be input either in P form or F form 
(2) D(s) can be input either in P form or F form 
(3) K values must be all positive or all negative. If 
both are desired, two separate runs must be made. Also, the 
maximum gain value cannot be zero. 
(4) An option card must be included to indicate whether or not 
a particular region of the root locus is to be drawn (zoom capa- 
bility). A blank option card implies no option isdesired. Note 
that selecting a specific region improves the accuracy of the plot. 
The last card tells the computer to plot only 


the roots locus in the rectangle in the s plane defined by: 


o . < Re[s] < 9 
- — "m 


min ax 
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En" - < 
Den t Im[s] <W 


max 


as illustrated in Figure 3-1. 





Standard Magnified 


Fig 3-1 Magnified and Standard Root Locus 


Thus a standard root locus plot is obtained by leaving the 
option card blank, while a magnified root locus is plotted 
by punching a "1" in the first column and specifying the 
minimum and maximum values of o and w. The input formats 
for RTLOC are given in Table Il. 
b. Output 

The problem identification is given, followed by 

the numerator and denominator polynomials, both in factored 


and 'ascending coefficients' form, and the minimum and maximum 
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BEE CA 


mpate Description 


Format 


Problem Identification 5A4 


letter P or F (for P 


nmm and F form), 
Order of N(s) 


Enters s) in format 
specified on previous 


card 


Letter Por F (for P 


(Sao, 


form and F form), 
Order of D(s) 


Enter D(s) 
specified on previous 


card 


Minimum value of gain, 
maximum value of gain 


(4 0) 


(< 10) 


in Format 


Al, I2 


SE10.0 


Al, I2 


8E1100 


HEP 


No option - blank card je, 9% 
Option # 0 


minimum 
maximum 
minimum 
maximum 


Table 


value 
value 
value 
value 


Ir =- Input Format Table for RTLOC 


Of ©, 
EE or 
GT 
Of w 


S ESOS 
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Columns Used 


I 


1% 2-3 


1-10, 11-20, 
etc. 


1-10, 11-20, 
etc. 


ege 


E Zoe 
201—307 pd 
41-50 





gains. The roots' real and imaginary parts are then listed 
as the gain varies from its minimum to maximum value. 
Finally the root locus plot is printed out. Note that the 
graph produced has square grids so that the true angles can 
be measured. 

This is normally a class B program and time = 2 
should be specified on the JOB card. 

c. Example 
Obtain the root locus of the following feedback 


control system: 






K(s + 1) 
Teen 






s2+ 4s+5 
3 2 
s” +7s +20s +50 






Fig 3-2 Feedback Control System for RTLOC Test. 


The equation for which the roots are to be found is then: 


(s+]) (s^ 4 4s 4 5) 
actus s3 + 7s* + 20s + 50) 


Il 
© 


ect KK 


DIS 





It agrees with the RTLOC structure so one can proceed further. 


N(s) (sal) (s +2441) (s + 2-741) 


and 


fede ede0s- + llos? + dle + 105° + 5° 


D(s) 


Here it is easier to enter N(s) in factored form and D(s) 
as an ascending polynomial. 

As a first quess, the range of variation of the 
gain is chosen to be from 0.0 to 100.0 and since the expected 
plot is unknown, no option is taken. 

This completes the work. The computer does the 
rest provided the cards are punched as follows: 

// (standard OS JOB card), TIME=2 
WWEXEC-LINCON 

WIbINK.SYSINADD,* 

NS TNCLUDE SYSLIB(RTLOC) 

/* 

//GO.SYSIN.DD.* 

EXBOC TEST 


F03 


0.0 0.0 150. Tio: 4l. 10. 
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0.0 100; 
(blank card) 
/* 

The results appear in Figs. 3-3A and 3-3B. Note 
that the user should mark the open-loop poles and zeroes for 
easier interpretation. 

2. Root Locus (PRTLOC) 

As the name indicates, this subprogram is a modified 
version of RTLOC. It calculates the roots of a polynomial 
and plots them. The method to input the data differs slightly 
but the ultimate goal remains the same. 

ARE DE 

This subprogram can only be used under Mode Two 
of operation. The coefficients of the polynomial must be 
entered using a simple subroutine called RPOL(C,G) which must 
be typed as follows: 

SUBROUTINE RPOL(C,G) 

DIMENSION C(20) 

C(1) = Énct (G) 


C2) = fnct (G) 


END 


where n = order of the eauation. 
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RUCT LOCUS PROGRAM 
PROSLEM IDENTIFICATION - RTLCC TEST 


RBSTSSHSS*ae PTs Seals CASS Ss SSE SEARS*SS SST SKEsERaAteaesesesvxnee 


NUMERATOR COEFFICIENTS IN ASCENDING POWERS OF S 


39 GAIN iz 7.126€ C1 
pCITS Se 
PEAL PL£2T IMAGO, FART 
-64533t 90 C.Q 
«.e25F-21 -3.Sol1E 2C 
4.^0527-0l1 3,901" CC 
2.2388 79 6.9 
-)1.0?9E 0) -7.53 1 E- IM 
-1.024E 92 7.9315F=91 
de GAIM « "S.E9CE 0? 
ROOTS Arte 
REAL PEAT [“46,. FAST 
5.57TEE-IL 4.1C8E GC 
-=6.655£ Jv Ce) 
—-2e21¢2— 00 4.0 
=} .l 92E 0) -7.277-91 
-1,103£ 039 7.2/7E-01 


Figure 3-3A Root Locus Test - Numerical Output 
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5.000 9.000 5.006 1.000 
CPEN-LOOP ZEROES 
REAL FART IMAG. PART 
-2.0C0E 09 -1.000E 00 
-2.000E QO 1.000€ 00 
-1.000€ 00 0.0 
DENCMINATOR COEFFICIENTS IN ASCENDING POWERS OF S 
0.0 0.0 150.000 110.000 41.000 10.000 1.000 
CPEI-LOOP POLES | 
EAL PART IMAGo PART 
-1.009E 00 -23.090E€ 00 
-1.000E 00 3.000€ 00 
-5.0U0€ 00 0.0 
3.000€ 00 0.0 
0.0 0.0 
0.0 0.0 
MIN. GAIN = 0.0 MAX. GAIN = 1.00E 02 
PESO LETTURE TTI TTT TIE 
1 GAIN = 0.0 
ROOTS APE 
REAL PART IM4Ge PART 
1.900€ 00 3.000€ 00 
-5.900€ 00 0.0 
0.0 0.0 
0.0 0.0 
2 GAIN = 5,.750E-02 
ROOTS APE 
REAL PART IMAG. PART 
-9.985E-0l1 -3.U00€ 00 
-9.985E-O01 3.0UCE 00 
-5.002€ 09 0.0 
-2.999EFE 09 0.0 
-1.023°-03 -4.37T7E-02 
-1.023E-03 4.377E€-02 
3 GAIN = 1.236E-Ol 
4 POOTS AFE 
REAL PART ING. PART 
-29.369t-01 -3.001* OQ 
-9.369E-O01 3.U01É 90 
-5.020*E 00 0.0 
= 2.1998 -03 -6-416E-02 
-2.199E£-03 6.+16ÉE-02 
D GAIN = 1.9Y7E-01 
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Hus... Cc (ntl) coefficients of the polynomial, in ascending 
order. Note that the coefficient of the 
highest order term C(n+1) must be 1.0 and 


need not be entered. 


mace (G) defining coefficient equation in terms of 


G, the gain. The function could very well 
be a constant only. 
The remaining data, i.e., the problem identifica- 


tion, range of gain values and option,are entered as follows: 


ENTRY Input Description Format Columns Used 

1 Problem identification, 544, 1 - 20 
Order of the polynomial, 12 DI 2 

2 minimum value of gain, EISE) jc oM 29 
maximum value of gain 

(4 0) 

3 {no option = blank card} Il, 9x, 1, 11-20, 21-30, 31-40 

option # 0, 8E10.0 41-50 


minimum value of o, 
maximum value of o, 
minimum value of w, 


maximum value of uw. 


Table III - Input Format Table for PRTLOC 


Here again, the gain values must be either all 


positive or all negative and the maximum gain cannot equal 


Zero. 
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identification must contain the polynomial order in columns 


E22. 


The first card, in addition to the usual problem 


The last card is used to 


not a portion only of the root locus 


plotted. 


If the option is selected, 


zero is punched in the first column, 


11-50 by the parameters definina the 
be blown up (see example). 


the card is left blank. 


indicate 
is to be 
a number 


followed 


rectangular portion to 


whether or 
refined and 
greater than 


in columns 


If this option is not desired, 


Note that this version permits us 


to find the roots of any characteristic equation with a single 


varying parameter G. 


b. 


the maximum gain are first listed out for future reference. 
Next, the root values are given as the gain varies and the 


meot locus plotted. 


Output 


The problem identification and the minimum and 


JOB card should be "time = 2". 


Co 


[3], one comes across the following characteristic equation 


Example 


The execution time to be included on the 


While trying to solve problem 7.26 in Shinners 


for part of the system: 


4 


S + 9.15s 


3 


+ (1.32 + 20K,)s“ + (26K, - .15)s + (6K, + 0.675) 


At this point the root locus is desired to determine what 


value of K 


2 


is required to satisfy some criterion. Since 
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the characteristic equation is specified explicitly, PRTLOC 
is selected. 

First the coefficients are sorted out and written 
as functions of G where G is equal to K,. 


s**0 coefficient : C(]l) 0.675 + 6.*G 


S**1 coefficient : C(2) -0.15 + 26.*G 


S**2 coefficient : C(3) 


1.32 + 20.*G 

S**3 coefficient: C(4) = 9.15 
Note that the coefficient of the highest order term is always 
taken as 1.0 and need not be included. The above data is to 
be entered by writing the subroutine RPOL(C,G). 

The order of the equation is 04. The range of 
gain values to be investigated is from 0.0 to 20.0 and since 
no refined plot is desired at this point the last card is 
a blank card. 

The following cards then constitute the entire 
deck to be input to the computer: 

// (standard OS JOB card) ,TIME=2 
77 EXEC LINCONF 
WIiGORT.SYSIN DD, * 
SUBROUTINE RPOL(C,G) 
DIMENSION C(20) 
C(1)=0.675+6.0*G 
C(2)==0:15+26.0*G 
ai =- 132720.: 7G 
C(4)=9.15 
RETURN 


END 
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Ss 

gratis v SS RN DD., 3 

an INCLUDE, SYSLIB (PRTLOC) 

n 1 ENTRYAPRTLOC 

/* 

//GO.SYSYN. DD. * 

PRTLOC TEST ONE 04 

0.0 20.0 

(blank card) 

A 

The results obtained with this first run as are shown in Figs. 
3-4A and 3-4B. However they do not permit us to evaluate the 
gain precisely enough and a second run is made, this time 
using the option. The rectangular portion where magnification 


1s desired is defined by: 


g = ze, 
min 
g = 5. 
max 
: = -]. 
“min 
= 5, 
“max 


Note that this option not only concerns the plotting but also 
produces a larger number of gain values. Thus, in order not 
to have too many values listed out unnecessarily, it is good 
practice to re-specify the range. It was decided to change 


metto vary from 0.0 to 10.0. 


47 





ROOT LUCUS PROGRAM 
PROBLE" TOEMTIFICATION = PRILCC TEST ONE 


*2$9295*9$9**t**t$xt*t$9?95$93»*t959*99595599à4*45kftt*55299 
MIM. GAIN » 0.0 MAX. GAIN = 2.00E 01 
SSSCSASSSEC SCS ec SAMe ees en CHS HFEVHSKVASH SAAR HERE 
1 GA IN = 0.C 


COTS Aft 
ART IMAG. PART 


GAIN = 5.750E-02 


PCOTS ARE 
REAL PART IMG, PART 


-8.888- Q9) 0.0 
1.072F-01 -4, Vi 7E-01 
1.072€-01 4.737E-01 

-9.168E-01 00 


3 Gäll = 1.4.236E€-01 
KOCITS ARE 
REAL FART ING. FART 
-5.755€ 99 0.0 
2.868: =J2 -5.972E-01 
2. ALBE-02  5.9T72E-01 
-4,521E-01 0.0 
s GàIM = l.4997E-01 
kQITS AQE 
REAL PA2T IMG. PART 
-8,597F 0) 0.0 
6.193E-92 -T.10UCE-91 
-6.193E-02  1.100F-01 
-4.28 37-01 0.0 
5 BAIN = 2.871-01 
ROOTS ARE 
PEAL PART IMAG. PART 
-8.41 0t 9) 0.0 
~1.66¢F-01 -8.201E-01 
-1.666€E—21 B.2ULE-Ol1 
+ =%.CI2£-01. 0.0 
6 GAIN = )3.877E-OL 
ROWS Alc 
REAL PART IMAGO, PART 
-8.194€ 0) 0.9 
-2.438d^-01 -9.278F-01 
-2.U8B8E-01] 9.-278F-Ol 
-23.586F-ul 0.0 
7 GAIN = 5,033E-01 


RCTS APE 
REAL PART IMaG, PART 


-7.919F 09 0.0 
~4.337¢-01 -1.03 
=4.331E-01 1.03 
-3.T23E-01L 0.0 


23 GAIN = 1.€3CE Ol 


PNOTS £P€ 
REAL PART  IMAG. PART 


-3.919E YY -1.73€F 0l 
sic 9iut 00 1.136€ Cl 


21.228€ CU 0.0 
-3.)26€-C1 0.9 
23 GAIN = LeESLE Ol 
FGOTS ERE 
REAL PART IMAG. FART 
-» 1.024€ 0U Del 
-3.412É 00 -1.375€ Ol 
-—3,912E 09 1.375€ 01 
-3.023€-01 0.0 


Figure 3-4A PRTLOC Test One - Numerical Output 
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Since the characteristic ecuation did not change, 
only the data deck is to be modified. These last three 
cards are given below. 
PRTLOC,TEST TWO 04 
0.0 10.0 
L x50 50 O 5.0 
This magnified portion of the root locus is 
presented in Fig. 3-5A and 3-5B 
3. Frequency Resvonse (FRESP) 
This subprogram determines the frequency response of 


a rational transfer function 


G(s) = K 





and plots the response in the form of a Bode or/and Nyquist 
diagram, as specified. 
a. Input 

The problem identification, the gain and the two 
polynomials N(s) and D(s) are entered followed by the minimum 
and the maximum radian frequency values, the number of fre- 
quencv values to be used (smaller or equal to 500), the inter- 
polation and discrete value options, the Bode plots and the 
Nyquist diagrams options and, only if required, the discrete 
frequency values. 

It might look complex at first, but the subprogram 
is very simple to use and the results obtained are auite good. 


The routine is accessed under Mode One or Mode Three. The 
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ROCT LOCUS PROGRAM 
FRUALE* IOENTIFICATIGN - FRTLCC TEST TWO 


99$9*9*99S*t9€4€994$9*2*t994*$9294€6959$4$459*$999»99&9895*9*$58$2924 
MIN. GAIN = 0.0 MAX, GAIN x 1.00€ Ol 
$€*$*€*X€69*€99gt$495y$9$5$9uv*394t*5$*9$99$99$*299$99$**$979t295»$$ 


OPTION HAS BEEN TAKEN 


SIGMA MIN = -5.,00F 00 SIGMA MAX = 5.00€ 00 
OMEGA MIN = -1.0CE 00 OMEGA MAX 2 5.00E 09 
$$99949999*992991$9929$$*2»29$3$*9$9999$$99$$92$*»9$92*9*9*$9$2$2829 

1 GAIN = 0.0 


ACOTS ARE 
REAL PART IMAG. PART 


-9.001E OO 0.0 

1. 7526-31 lee sì 
~4.99TE=-01 0.0 

2 GAIN = 2.080E-02 


ROOTS 
REAL PART 


-8.960E 09 


AR 

IMAG. PART 

0 
1.507=-=01 = 

0 


1.9507 -91 
-4.913E-91 


T 
J 
l 
L 
1 

AIN = 6.493E-02 

TS AFE 
REEL r 
J 
e 
2 
L 


-8.87 
9.303c 


-4.1% 
5 GAIN =  8.833E-02 


IMAGO, PART 


Case tal 


PCOTS 
REAL PART 
-5.825F 09 


6 GAIN = L.LI27E-OLl 


OITS ARE 
A I4AG, PART 


To GAIN = 9.3318 CC 


FOTS ARE 


KCAL PART [MAG. PART 
anna) =-1,.3C2E II 
-1.048E NGO 0.0 
ze A eh E sei L 0-9 


Figure 3-5A  PRTLOC Test Two - Numerical Output 
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input format table and the example that follows demonstrate 


the procedure to be used. 


Entry 
1l 


2 


8 
(if and 
only if 
BPLIONT 
= 001 


Input Description 


Problem identification 


The gain K 


letter P or F (for P form 


Or F form), 
N(s) < 10 


Enter N(s) 


order of 


in format 


specified on the previous 


card 


letter P or F (for P form 


OPEP form), 
D(s) « 10 


enter D(s) 


order of 


in format 


specified on the 
previous card 


minimum radian frequency (#0), 


maximum radian frequency, 
number of frequency values 


to be used 


eptien I: 


(<500) 


logarithmic 


interpolation 


= 000 


discrete values 


supplied 


= 001 


linear interpolation = 


Operon B: 


option N: 


Bode plot = 
no Bođe plot 


Nyquist plot 


no Nyquist plot 


= 001 


002 


000 


discrete frequency values 


Table IV - Input Format Table for FRESP 
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001 


000 


Format 


5A4 


8E10.0 


Al,I2 


8El0.0 


Al, I2 


8E10.0 


2E10.0 


73 


I3 


ES 


8E10.0 


Columns Used 
1-10 
1-10 


E o 


l-10, 
etc. 


11-20, 


e 


1-10 , 
etc. 


117220, 


LALO; 
11-20, 
21-23; 


24-26, 


21229 


20232 


1-10, 
etc. 


ep 





Here the option card is a bit complex, but it 
provides great flexibility. The following ideas should help 
clarify the concept: 

(1) The first three entries specify the range and the 
number of data points for the Bode and/or Nyquist plot. One 
must recall that the Bode magnitude plot is log-log; the Bode 
phase plot is log-linear (angles in degree) while the Nyquist 
is a polar plot. Thus, minimum and maximum radian frequency 
values should be carefully chosen. For example, w_._ = 0.01 


man 


and e 100 could be a good choice in a given problem 
while being absurd for another one. 
(2) Option I specifies the type of interpolation to be used 
to generate the values between the minimum and maximum frequency. 
If Option I = 000, logarithmic interpolation is used to 
select the frequency value. Either plot 


can be obtained while specifying this option. 


00l, the user must enter on the following 


IEF Option I 
cards the frequency values for which 
he wants G(jw) to be evaluated. The 
number of frequency values must 
again be less or equal to 500. No 
plot can be obtained when this option 
is selected, only tabular outputs. 

If Option I = 002, linear interpolation is used to select 
the frequency values for which G(jw) 
is to be computed. Only the Nyquist 
plot can be obtained when this option 


is used. 
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(3) Option B indicates whether or not Bode diagrams are 
to be drawn. 


If Option B = 000, Bode plots will be output 


If Option B = 001, Bode plots will not be output 

(4) Similarly, option N is used to specify whether a 

Nyquist plot is desired or not. 
If option N = 000, it is desired 
If option N = 001, it is not desired. 

(5) The options card is not followed by any card except 
when option I is equal to 001. If this is the case, the 
frequency values must be entered using an 8E10.0 format. 
Note that an option card containing only the minimum and the 
maximum frequency values and the number of points to be 
evaluated indicates that both Bode and Nyquist plot are 
desired. 

D. tout 

The problem identification, the value of the gain, 
the coefficients of the polynomials N(s) and D(s) as well as 
their roots are listed for reference. Next, the radian fre- 
quency, the real and imaginary part of G(jw), the magnitude 
IG(jw) |, the magnitude in db, the phase in radians and the 
phase in degrees are printed out in tabular form for the 
indicated number of frequency values (smaller or equal to 
500). 

If option B = 000 has been selected, the magnitude 


and phase Bode diagrams are given. Note that the phase angles 
are always normalized to lie between -180 and +180°. 
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If option N has been requested, the Nyquist 
Braecram 1s plotted with the points linearly, or logarithmically, 
spaced out. 

Normally the CPU time required to run the program 
is less than 20 seconds (class A). 


c. Example One 


The Bode plot for the loop transfer function is to 


be obtained for the following systen: 





60 (s +50) 
s(s +10) 


Fig. 3-6 Compensated Control System for FRESP Test 


The first step is to define G(s). Here it is simply 


(s + 50) 


Be $ 7 S0lStS¥I0) (s¥ 20)! 





The gain is 60 and since both N(s) and D(s) are already 
factored, they can best be entered using F form. The minimum 
and maximum frequencies are arbitrarily chosen to be 0.1 and 
100, respectively. The number of frequency values for which 
G(jw) is to be evaluated is 50. Since a Bode diagram is 


desired, option I must equal 000 (logarithmic interpolation) 
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and option B also equals 000. In this case a Nyquist plot 
is not desired and option N is entered as 001. 
The control cards and data deck to run the sub- 
program are then: 
// (standard OS JOB card) 
// »EXEC „LINCON 
VABTINK.SYSIN,DD,* 
NINE LUDE SYSLIB(FRESP) 
p^ 
WiGO.SYSIN DD * 
FRESP TEST ONE 
60. 
FOL 
50. 


E03 


0.1 100. 050000000001 
SE 

The Bode diagrams are shown in Fig. 3-7 (A-C). Note 
that the phase versus frequency diagram presented is not in 
error but simply due to the fact that the angle values are 
normalized to be within -180 and +180 degrees. This is 
useful since it permits one to rapidly determine where the 


-180° crossing occurs. 
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PEASE (NEG) 


(RADI 


PHASE 


IN DA 


MAG 


MAGNITUDE 


000C00€ 00 


FCWERS OF S 


POWERS OF S 


IMAGINARY PART 


IN ASCENDING 


2.600000€ 02 3.000000F Ol 1 


ENUMIMATOR ROOTS ARE 


In ASCFNDING 


PART 
PART 


INLG, 


IMAG. 


SE 
CETLCN - FRESP TEST ONE 
REAL PARI 


yy 


LF 
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PF 


COEFFICIENTS 


Ul 0.0 


DEMANA TOR 
Us 


- 
= 


ke 
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o« C90000E C1 
Pt AL. PART 
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d. Example Two 
The problem to be solved now requires that only 
a Nyquist plot be obtained for the following open-loop 


transfer function of a compensated system: 


E. sd) ls +50 1320 — 
(s +7) (s +0.015) (s +1) (s+2)s 


The frequency range is selected to be from 0.2 to 10.0 and 
calculations are to be carried out for twenty-five frequency 
values using logarithmic interpolation. The:computer deck is 
then: 

// (standard OS JOB card) 

Aa EXEC, LINCON 

EEXSINK.SYSINADD,* 

RRENCLUDE.SYSLIB(FRESP) 

ye 

//GO.SYSIN,DD,* 

FRESP TEST TWO 

20.0 

F02 

E. 15 

nz 7 


F05 


6l 





p 


0.2 10.0 025000001000 
/* 
Results are shown in Figs. 3-8A and 3-8B. 
4. Partial Fraction Expansion (PRFEXP) 
This subprogram performs the partial fraction 


expansion of the ratio of two polynomials of the form 


D(s) ' (degree of N(s) < 
degree of D(s)) 
Elo Aba 

The problem identification, the gain value K and 


the polynomials N(s) and D(s) are entered according to the 


following input format table: 


ENTRY Input Description Format Columns Used 

l Problem identification 5A4 1-20 

2 gain value K 8F10.3 WETO 

3 letter P or F (for P form Al, 12 1, 2-3 
and = form), order of 
N(s) < 10 

4 enter N(s) in form specified 8r10.0 1-10, 10-11, 
on previous card etc. 

5 Torre Por F (for P form Al, I2 1, 2-3 
and F form), order of 
Dis) < 10 

6 enter D(s) in form specified 8F10.0 l-10, 10-11, 
on previous card etc. 


Table V - Input Format Table for PRFEXP 
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Figure 3-8A FRESP Test Two - Tabular Output 
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D(s) must not have multiple complex roots for the subprogram 
to work. If it does, a message is printed and the problem 
terminates at that point. Note that D(s) may have multiple 
real roots though. 
DP OUEPUt 
The problem identification and the gain value are 
listed followed by the numerator and the denominator in both 
factored and unfactored forms. For the denominator, 
each root value is listed once only with its multiplicity 
indicated. Note that roots are considered equal if their real 
and imaginary parts do not differ by more than 0.005. The 
example presented in c, illustrates how to deal with multi- 
plicity of roots in the interpretation of the results. The 
residue matrix real and imaginary parts NER hen given. 
This subprogram can be run as a class A job. 
c. Example 
The partial fraction expansions of the following 
rational functions are to be performed: 
20 


(a) 


656 4 752s 4 264s? 4 28s? 4 s^ 


2+s 
m wir 20085 


2+4s+3s +85 


N(s) and D(s) are entered using both the F and the P forms and 
the partial fraction expansions of the two polynomial ratios can 


easily be obtained in a single run by stacking the data deck. 
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The computer cards are: 
// (standard OS JOB card) 
// EXEC LINCON 
gees VG CN DDx* 
INCLUDE  SYSLIB (PRFEXP) 
/* 
//GO.SYSIN,DD,* 
PARTIAL FRACTION A 
20. 
POO 
ie 0 
P04 
056. 752: 264. 20 TE 


PARTIAL FRACTION B 


/* 
and the solutions are presented in Fias. 3-9A and 3-9B. 


Interpretation of these results gives: 


| _ .0139 4 j.0124 , .0139 - j.0124 
Gye Partial fraction A = rar AT + S 1231-47 


-.0278 n .1667 
s+2 (s+2) * 
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PARTIAL FRACTION son 
PRCSLEM IDENTIFICATION = PARTIAL FRACTION A 


NUMFRATGR GAIN s 2.0COE C1 
NUMERATCR CCcFF. - IN ASCENCINC PQXERS 

1.029€ 09 
DENOMINATSR CCEFF. - IN ASCEACIAG POWERS 

6.360€ 02 7.520€ 02 2.540€ 02 2.8008 O1 1.00CE 00 
SRVSEKTE SAK TH SKAKH VIR SK SS SAAB EKE ADEE HEAER 


CENCHINATCR RCOTS 


REL PART IMZ Çe PART MULTIPLICITY 
-1.15999925 01-4.4721594é 00 1 
-l.19$$9$92€ Ol 4.47215943 00 1 
-2.0011177¿ 00 0.0 2 


*59:4ux5w9415u€ $535 2 41 36 09 OA c e ade de e a a a Gov Go ad A 


RESIDUI MATRIX = REAL PART 


1734392205*-02 
0 3892505¢-02 
=2.1783011=-02 1.6€£5752F-01 


RESICUS MATRIX -— IMAG. PART 
1.2423638- -22 
=1.24236338-02 

0.0 9.0 


BARACCA AAA AS A A A ae a B 


Figure 3-9A Partial Fraction Expansion A 


PARTIAL FRACTICN EXPANSION 
PROSLE4 IDENTIFICATICN - PARTIAL FRACTICN B 
AUMERATOR GAIN = 1.000€ 00 


NUMERATOR COEFF., — IN ASCENDING POWERS 
2.930€ 07 1.0932 Cu 
MJMEKATOR RCOTS 


REAL PART IMAG. FART 
-2 GJ000G0F 00 0.0 


DENSMINATOR COEFF. - IN ASCENOING POWERS 
2.900€ 00 4-000F 00 3.00CE 00 1.000E 00 
SBSFSRAVISSK MAA 185 39 5 308 9 9 49 0 X90 23 Rd NY 
CENCMINATOR RGOCTS 
REAL PART IM 


A .T_ FULTIPLICITY 
71.079205 90-1.29 1 


n 
=1 259000632 CC LS G 1 
29 :95999543-01 0.0 1 


BEREREKABUEREANTZLKGTEREZEREREHETEIDIKER HH TH AH a A ad 
RESIOUE MATRIX = REAL PART 
-4.,99999C5*2-01. 
—. 95595592057 -01 

9.995994CIE-01 
RESIDUE MATRIX - IMAG. PART 
—.9999923£5-01 

Q.0 

*30935X 99150): 915 39 € Xp 03 00300098 o a e a ad a OG X 


Figure 3-9B Partial Fraction Expansion B 
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Note that the second residue appearing in the output belongs 
2 ane 
to (s+2) . If a multiplicity three had been the case, a 


third residue would have been the numerator of a cubic. 


; LS SS 1 
2 Partia Ce ee eps <+1 
(2) lE ractron B -stltj * sH-j * srl 


TROOS Of a Polynomial (Roots) 


This subprogram finds the roots of a polynomial of 
degree less or equal to twenty. 

SE Input 

The first data card contains the problem SE 

fication in the first twenty columns and the polynomial order 
in columns 21-22 (format I2). On the next card(s) the poly- 
nomial coefficients starting with the lowest order term are 
entered (format 8E10.0). These two entries are repeated for 
every polynomial to be factored. Note that the highest order 


term coefficient must be unity. 


Entry Input Description Format Columns Used 
1 Problem identification, 5A4 1-20 
polynomial order 12 21-22 
2 polynomial coefficients SEA 1-10, 11-20, 
in ascending order ENER 


(highest order term 
coefficient being One) 
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ba Output 
The problem identification and the polynomial 
coefficients are listed for reference. The roots real and 
imaginary parts are then printed. 
c. Example 


The following polynomials are to be factored: 


The computer cards are then: 
//(standard OS JOB card) 

/ / ABXEC ALINCON 
estate SS VG TN DÉI. $ 

na a INCLUDE SYSLIB (ROOTS) 

/* 

EGO. SYSINADD, * 

Roots test one 03 

de. O 0.0 0.0 1:0 

Roots test two 04 

1 -5.0 E20 1.0 O 
Roots test three 

0.0 1.0 1.0 0.0 050 1.0 


* 
The result is shown in Figure 3-10. 


C. TIME RESPONSE AND MATRIX MANIPULATION SUBPROGRAMS 


These three subprograms permit a user to analyze linear 


control systems for rational and graphical time response and 
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also provide matrix manipulation to easily solve for deter- 
minants, inverses, state transition and resolvent matrices, 
eigenvalues and characteristic polynomials. 

The control system must be linear and represented in state 


variable form as [1] 


x(t)=A x(t) + b ult) 

_ T 
u(t)=K[r(t) - k"x(t)] 
y (t)=c x(t) 


where u(t), r(t) and K are scalar and the system order is less 


or equal to ten. In block diagram form, the matrix system can 


be represented as 


y(t) 







controller 
K 





Fig 3-11 Linear Control System Block Diagram. 


The diagram readily shows all the elements to be provided for 
the study of any given system. For instance one can view that 
setting k* = 0 gives an open-loop system and that unforced sys- 


tem analysis can be done by simply letting r(t) = 0. 
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1. Basic Matrix Manipulation (BASMAT) 


This subprogram is used to perform various calculations 
associated with the plant matrix A of a given linear control 
system. It is a class A job and must be run under Mode One 
or Mode Three. 

a. Input 

The problem identification and the dimension of A 
are given on the first card. Next the A matrix is entered, 
one row at a time using an 8E10.5 format. Thus, if the dimen- 
sion of the matrix is eight or less, one row per card. Other- 
wise the 9th and/or 10th elements appear on a second card and 
the rule becomes one row per two cards. The last card indi- 
cates what matrix operations are to be performed. The key to 


obtain the proper results is explained after the input format 


table. 
Entry Input Description Format Columns Used 
1 Problem identification 5A4, I2 l, 2-3 
dimension of A(N < 10) 
2 A(N xN) matrix 8E10.5 l-10,11-20,etc. 
(one row per card 
for N < 8; one row per two 8E10.5 1-105259 -20,etc. 
cards for N»8) 
3 option det 0, determinant Il B 
desired 
1, determinant Il 


not desired 


option inv inverse desired Il PAP 


— © 
= 


, inverse not 
desired 


72 





Entry Input Description Format Columns Used 


option phi (s) + 0, $ (s) desired Ib 3 


L; (s) net 
i desired 


epedon C.E. 0, characteristic Il 4, 
polynomial 
desired 


l, characteristic 
polynomial 
not desired 


option eigen 0, eigenvalues Il oe 
desired 


l, eigenvalues 
not desired 


option phi (t) 0, d(t) desired Il 6. 
l0 E)  HOt 
P desired 


Table VI - Input Format Table for BASMAT 


Thus, a zero indicates that the computation is desired while 
a number from l to 9 informs that the listed operation is not 
to be performed. Six zeros or a blank card would result in an 
output that contains the A matrix determinant, inverse, resolvent, 
characteristic polynomial, eigenvalues and state transition matrix. 
be Output 
The problem identification and the A matrix are 

listed first. Then the result of each operation selected on 
the option card is printed as follows: 

(1) det (A) - a scalar 

(2) a> - a matrix presented one row at a time 

(3), (4) resolvent matrix and characteristic polynomial. 

ls (s) = [si -A]7?, which is called the resolvent matrix, 


is the Laplace" transform of the state transition matrix 


p(t) = e2, 
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une coefficient matrix of the numerator of the 
resolvent matrix appears first, followed by the character- 
istic polynomial in ascending powers of s. 

(5) Eigenvalues - listed indicating the real and imaginary 
parts 

(6) Time domain state transition matrix - $(t) (see part c, 
example two). 

The subprogram is restricted by the fact that $ (t) 
cannot be calculated if eigenvalues are multiple. If a situa- 
tion where the state transition matrix is requested where 
eigenvalues are not simple, a message is printed (see part 
c, example one) and the computer goes to the next problem. 
Note that eigenvalues are considered to be identical if their 
real parts and their imaginary parts differ by less than 
0.005. 

c. Examples 
(1) Example One 
The resolvent matrix (SI - Si and the state 


transition matrix $(t) are to be found for the plant matrix 


0 d 0 
A = 0 0 Il 
0 0 =2 


Thus, the matrix has dimension N = 3. The options are set 


equal to the following values: 
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option det = 1 determinant value is not desired 


option inv = 1 inverse n° is not to be calculated 
option phi(s) = 0 $ (s) is desired 

option C.E. = 0 characteristic polynomial is desired 
option eigen = 0 eigenvalues are to be computed 
option phi(t) = 0 $(t) is desired. 


The computer card deck is then: 
// (standard OS JOB card) 

// ,EXEC, LINCON 
//LINK.SYSIN_DD_* 

_ INCLUDE, SYSLIB (BASMAT) 

FR 

27/G0.5SYSIN, DD _* 


BASMAT TEST ONE 03 


0.0 1.0 0.0 
0.0 0.0 ie 
0.0 0.0 -2.0 
110000 
/* 
The computer results shown in Fig 3-12 can be interpreted as 
follows: 
642s s+2 È E E ———- 
S S +2) 
1 > S s (s 
$(s) = EX v 0 s +2s ss = 0 I l 
~ (e +2s°) 2 S s (s+2) 
0 i 0 S 
l 
p 0 sr? 
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ZASIC as en : 
PROBLEM IUEKTIFICATION- BASMAT TEST CNE 


BU LFEERAEE KA TRERPFERETBEBERELTEHTCEHP PETE RL FR RE 


THE A MATRIX 


@.JUCQQIVE 00 è e 5 CNECULE 00 1,0C00000E£ 00" 

0.0 vos CÉCCLIOCE OO 1.0000000€ 00 

0.) ~4,UCCCCCCE OO 2«0000000F 00 
$$$ 6e»2944«v$$4$&à5$9*55a8$9$»v*nu$»28522952519*8a999$v0*5à€**t*599 


THE MATPIX CCEFFICIENTS CF THE NUMFRATNR CF THE PESOLVENT MATRIX 
THE MATRIX COEFFICIENT CF S442 


de case 00 0.0 deo 
l.CCCCCCCE GO 0.0 
0: 5 0.0 1.000000NE 00 

THE MATRIX COSFFICIENT CF Ses) 

-8,9GC9090E OO ce eOCODCCCE 1,0009009F 00 
Qed -5«bLCoC.o0cE 0) Y,0C02020* 09 
"pue =6 JSOCCLULLE -1.0C00000° Ol 

TRF MATRIX CCEFFICIENT CF S**0 
1.59995 945c Ol -3.0CICCOCE 99 -4.0000000F 00 
GeV T.99SSERTE 0) -4.0000099€ 00 
0.0 I.6CCCCCOE Ol 2,3599985" 01 

titze roseto 


DBRBRERESZ RER PER EV EP RR HR I RE EFT 9249 


THE CHEFACTEFISTIC PCLYROMIAL — IN ASCENDING POWERS MF $ 


-6.39S99955E Q1 4.79%SS65E Ol -1.2000000E O1 1.0000000E 00 
SSSVSASAKTeKacrwewe eee st asst saseseeeesexaetew*eeuake 
TRE ELGENVALUFS OF Tre f MATRIX 
AEAL FAKT IMAGINERY PART 

€.00C1L345F N0 0.0 

@eCNCLID+WE 00 0.0 

3.9994(292t 09 0.0 
weewAElINGeee 


EIGZY.VALUCS MUST BE SIMPLE 
CALCULATIONS CANNGT BE COMPLETED FOR THIS PROBLEM 


9$9**:4$$*522*5299$259499$929957179?*4845*5$5*9599225959$5»* 


Figure 3-12  BASMAT Test One 
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The state transition matrix $ (t) cannot be obtained since the 
eigenvalues are not simple. 
(2) Example Two 
This second example shows the complete 
solution, i.e., determinant, inverse, resolvent matrix, 
characteristic polynomial, eigenvalues and state transition 


matrix, for a case where 


Since all the calculations are requested, the option card is 
left blank. The card deck is 

// (standard OS JOB card) 

// EXEC LINCON 

WBINK.SYSIN DD * 

__ INCLUDE SYSLIB(BASMAT) 

/* 

7/GO.SYSIN DD * 


BASMAT TEST TWO 03 


2.0 232 SE 
DE 3.4 pe 
E? HL IL 


(blank card) 


Ken 


n 








Results appear in Fig 3-13. Interpretation of these results 
is fairly straightforward. For instance, the first term of 


the resolvent matrix, ¢(s), is 


~~ 


BE - 4.5s - 4.07 


$14 (s) 
S - 6.55 - 8.54 s + 0.049 


and, similarly, the first term of the transition matrix, 
$(t), is 


+0. E 
$4, (6) -  (0.475e "0: 0097€ 4 q, 5595 1-13t 4 9,2960: 62, 


2. Rational Time Response (RTRESP) 
This subprogram may be used whenever it is desired to 
obtain the time response in closed form [1] of a linear control 


system described by the following set of equations: 


x(t) - A x(t) * b u(t) 

T 
eye = ONU) ko x6] 
y(t) = c x(t) 


The system can have any initial conditions x(t.) but the scalar 
forcing function r(t) must have a rational Laplace transform 


such that 


[r(t)] = R(s) = G = , where G is a constant, 
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Anal MATRIX bel uis 4 
Raat EM IVENTIFICATION- PASMAT TEST 


DERTDERDBERTATBERKRTUAREDBRTRELTELA IF PL EAT KERLE RR 


THE A YATE IX 


2.00260900* QO 2 1599054“ 09 

5,3444459 00 3.35555S6E 00 

8.93999998F-0Ol 1.059595$4€E OO 
THE OETEPMIKANT CF TFE FATRIX 


-4,9031435E-0? 
TRE INVERSE CFE THE MATRIX 


6.32964} 829E Ol =6.7341212€ 00 
~1.5915466e Ol 1.02041C55 00 
~5.20412249E Cl 4.485€157E 00 

FZRSILEZLIEPRL 


LE aa ieri rt 


THE MATRIX COEFFICIENTS CF TF 

THe MATRIX COFFFICISNT CF S**2 
1.000J990E 00 “3 .LSI17S78E-05 
del eCCocCCCc 00 

THE MATRIA CCEFFICIENT CF Sas] 

— @,495+59)F 00 2 61533553" 00 
5.059946594É 00 ~3.USISSS4E 09 
9.9999 933f-U1 1,C535554£ 00 

THE MATRIX COEFFICIENT CF Saa 

-4.0693921EF VO 3.2656652F-01 
7.7339733c-Cl 4,9334 4639-02 
2.549956%E 00 -2.19595231E-01 


BERURULN LET EKUB UAR ZUFUHR LTE FT TI Pd pria 

THe CHAFACTERISTIC PCLYMGMIAL - IN ASCEN 
4.9001485F-02 -8.5395904E 009 

+28 count tz tda vat LEUTE RB PAR RED IH BR U KL HT U I IF A RR 


THE EE cS OF THE A MATRIX 


QE ae IMAGINARY PART 
5.71350%2f.-03 JeU 
-1.1«50170* Ou Ved 
1.019359297t CO Jeu 
DAD NAAA ADUANA AAA RE TR I AAA rs 
THe ELEMENTS CF TRE STATE TRANSITION MAT 
THE MATRIX COEFFICIENT CF EXPL 5.713064F 
dr e 71545737=-01 EE 
-9.395217150*-02 71.693159 77-03 
-2.9662031E-01 2.94505€5St-02 
THE MATATX CCEFFICIENT CF ExPt-1.1256517E 
2.28€525^*-UCl -2.1693249F-01 
~5.0137544:.-0l 4.75042746E-01 
145933990t-01 -1.517371380C-01l 
THE MATRIX CCEFFICISNT CF EXPE 7.6 
2.5539075€-01 2,56101825-01 
- 509530637E-C1 $.164€215*-01 
1.4125095€-01 14225€515c-61 


TWO 


2.5C00000F 00 
7.0999994E 00 
1.0999994E 00 


ON 
t o e 

ooe 
AAV 
Q-Ou 
ANO 
DO «y 
NES 
C=O 
nnn 
ooo 
ume (NJ 


21:5258789E-05 
=1.32531896-05 
1,0C09000E 00 
2.5000N0CE 00 
"7.0999994E 00 
-5,3599996E 00 


DING POWERS OF 
-6.4999990E 00 


09)T 


199015 0037 


Figure 3-13 BASMAT Test Two 


po 


E NUMERATOR OF THE RESOLVENT MATRIX 


1. 0000000E 00 





and 


Sta e? e t 


N(s) = a 1 > s 


+ a 
0 | 


b +s + b. some... + SD 


eS) o ET 2 


milena m > £ > 0. 
Arrange the polynomials so the coefficients in the highest 
order terms of both N(s) and D(s) are unity and select the input 
gain G as required. 
In addition to the above, it is necessary that the 
total order of the system, i.e. order of D(s) plus dimension 
of A be smaller than or equal to ten. This limitation is not 
overly restrictive but must be taken into Ton: when handling 
large order systens. 
lr O 

The system matrices, feedback coefficients and the 
controller gain are entered immediately after the problem 
identification and system order card. The A matrix elements 
are presented one row at a time. The transpose control vector 
P the output vector c, the feedback coefficients ki, Ko, 
Ka, Ir k and the controller gain K are given using an 
8F10.4 format. 

Next the initial conditions x, (0), x5 (0), ti 
x, (0), the input gain G and the numerator and the denominator 
input polynomials are entered. Both N(s) and D(s) may be 


entered in factored (F) form or unfactored (P) form and it is 
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noted that the degree of D(s) must be strictly larger than the 
degree of N(s). 

It is suggested that a signal flow graph, or at 
least a matrix block diagram, be sketched before an attempt 
is made to run this subprogram. It does not take long to do 
so and much can be gained. 

The execution time for the subprogram is less than 


20 seconds for most cases (class A). 


Entry Input Description Format Columns Used 

D Problem identification, order 5A4, I2 1-20, 21-22 
of the system (N « 10) 

2 plant matrix A (one row per 8F10.4 1-10, 11-20, 
card if N < 8; one row per etc. 
two cards if N > 8) 

3 Control matrix EX x N) 8F10.4 1-10, 11-20, 
(on one card if N < 8; etc. 


two cards if N > 8) 


4 Output vector c(1l xN) 8F10.4 1-10, 11-20, 
(on one card if N « 8, etc. 
on two cards if N » 8) 
5 feedback coefficients k,,k 8F10.4 1-10, 11-20, 
Parsi eer, 
eee (on One card if etc. 
N < 8; on two cards if N > 8) 
6 Controller gain K 8F10.4 1-10 
7 Initial condition x, (0), 8F10.4 1-10, 11-20, 
«eni. X (0) (Gn one etc. 
card if N < 8, on two 
cards if N > 8) 
8 Input gain G 8F10.4 l-10 
9 Letter P or F (for P form Al, I2 1, 2-3 


or F form), polynomial 
order 2 < M 


81 





Entry 


10 


EL 


12 


for reference. 


both in factored and unfactored forms. 


Input Description 


Enter N(s) in format 
specified on the previous 
card. 


Letter P or F (for P form 
or F form), polynomial 
order M< 10 


Enter D(s) in format 
specified on the 
previous card. 


Table VII - Input Format Table for RTRESP 


b. Output Format 


Format 


8F10.4 


Al, I2 


8F10.4 


Columns Used 


l-10, 11-20, 
etc. 


lee 


1-10, 11-20, 
etc. 


All the information given as input is repeated 


The polynomials N(s) and D(s) are presented 


The rational time response of each component of 


the state vector x(t) and the scalar output y(t) are printed 


in pseudo-matrix form. 


elarify 


complex 


THE 


THE 


THE 


THE 


Here again a hypothetical example can 


the presentation. For a two-state problem, assuming 


poles and a step input, the computer output would look like: 


TIME RESPONSE OF THE STATE X(t) 


VECTOR COEFFICIENT OF EXP(A)T * COS(B)T 


X 


Zn 12 


VECTOR COEFFICIENT OF EXP (A)T * SIN(B)T 


SE Saa 


VECTOR COEFFICIENT OF EXP(0.0)T 


290 32 
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where X]}' X19! X51' X42! X31 and X35 are numbers. The result 


would be interpreted as: 


x, (t) X11 * exp(at) * cos(bt) -* Xo, * explat) * sin(bt) + x 


ET 


X Eexpiat) * cos(bt) + X 


12 * exp(at) * sin(bt) + x 


x, (t) 


22 32 


The procedure to obtain y(t) is the same. Note 
that if more than one output y(t) is desired, the subprogram 
must be rerun changing the c matrix each time. 

c. Example 
The open-loop rational time response is desired 


for 


K 
DI 





e 
EE 


A 
D 


The first step is to get the signal flow graph 


and state equations. 





Fig 3-14 Control System for RTRESP Test 
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me = x. (€) 


Sa 9 = -3.846x,(t)- 2.346x,(t) + .1923u(t) 
Wate) = r(t) 
y(t) = x4 (t) 


The data from the system is then: 


DNO ae 
A = 
È -3.846 -2.346 
pb? = [0.0 .1923] 
a [1.0 0.0] 
k = [0.0 0.0] 
K = 1.0 
x(0) = 0 


The system time response in closed form is required 


for a step input of magnitude 2. Thus 


R(s) = 


VIN 
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and the data are 


input gain = 2 
N(s) = 1 
D(s) = s 


The control and data cards to run the program are as follow: 
// (standard OS JOB card) 
// EXEC, LINCON 
//LINK.SYSIN,DD * 

~na INCLUDE, SYSLIB (RTRESP) 
ES 

//GO.SYSIN,DD,* 

RTRESP TEST 02 
0.0 1.0 

-3.836 -2.346 

0.0 0.923 


1.0 0.0 


FOL 
0.0 
/* 
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The results shown in Fig. 3-15 are interpreted as: 


x(t) NELLO I exp(-I-173t) * cos(1.57t) - .075 * exp rjr) 


sesinil 576) + 0.1 
x, (t) ANE d exp (-1.173t) * sin(1.57t) 


3. Graphical Time Response (GTRESP) 


The subprogram is a slightly modified version of the 
one presented by Melsa and Jones [1]. It still determines the 


time response of the closed loop system 


x(t) = A x(t) + b u(t) 


a(t) = Kir(t) - k” x(t)] 


y(t) C x (t) 


with initial conditions x (tà) and displays the results both 
in tabular and graphical forms. However, instead of having 
all the desired plots drawn on one graph only, it also produces 
one graph for every selected variable. 
The subprogram solves linear systems. It is a Class 
B job when graphical output is requested but reduces to a 
Class A job when tabular output only is to be listed. The 
subprogram must be accessed under Mode Two and requires an 


exterior subroutine to define the scalar forcing input r(t). 
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- VER HELLE I ALL UHR HA RI EHE 


RATICNAL TINE QESPCASE 
PROGLEM IUEATIFICATICN — RTRESP TEST 
sarna 


*s2$6$$92999299*9*555»5939 


TE A MAT2]x 


0.0 1.9000000 
-3.8£453997 -2.3459997 

THE AR MATRIX 

0.0 C.1923000 
THE C MATRIX 

1,3997300 C.0 
FEFOBACK CCEFF, 

3.0 0.0 


GAIN = 1.CCCINODE 00 
SHHSTHS SESS HHH EASA EZCKEKEKKKA TEE ESA EA CU 
INITTaL CNAMOITICNS - xX(Q) 
0.0 0.0 
RCAIN = 2.00G100CF 09 


NUMERATER PCLYNC"IàL OF R(S) - ASCENDING PCWFRS OF S 
1.C6222C60 ‘ 


EE RAR OF Ra z NING POWERS OF S 


ASCEN 
0090000 
DENSMISATIZ ROOTS ARE 

REAL FART IVAG. PART 

J.0 0.0 


GEEHRT I ITETLEEERIL TER I BETT RER OH IF HRATFT AT FE R 


TRE TIME RESOGUSE CF TRE STATF x(T) 


THE VECTOS COFFFICISNRT CF EXP(-1.172999C 00)T*COSI1.571647E 0017 
-9.16239226Q Q,2020201 
THE VECT9II CE=SFFICIENT CF EXP(-1.172999€£ OOJT*STIN(I1.571647T€ 00)T 
THE VECTCI COESFFICIENT CF ExF( 9.9 IT 
3.059936 -C.2000251 


FERE HT EHER TER TE THE TE TER LE DLR rete 


Une IE 22S5°CNSE CF TFF CUTPUT YITI 
ThE COM E320 EXP (-1.172999€ 00) T*CP^S(1.571647E 00)" 


-U. 10 y 
TRE a E CG EXP (-1.,172999E ONITSSTN(1.,571647F 90) T 
TRE CCEFFICIE NT DE EXPLI C.C IT 

0.6739 


SESSSCSHSH SSCS SSCHSSHRAS ESSE AIPKEKSEH SHH VES SISSHAeess 


Figure 3-15 RTRESP Test - Computer Output 
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œ Input 
The first element to be input is the forcing func- 
tion r(t). A short defining subroutine must be written in 
the following manner: 
SUBROUTINE RFIND (T,R) 
(FORTRAN statements defining r(t)) 
(Example: R = 2.5*T+SIN(4.2*T) 
RETURN 
END 
Next the remaining parameters are entered as a data deck which 
closely resembles the one for RTRESP. The problem identifica- 
tion and system order (N < 10) are given on the first card. 
Then the NxN plant matrix A, the single row matrix o the 
output matrix c, the feedback coefficient matrix kT, the con- 
troller gain K and the initial conditions x (t9) are presented 
as indicated on the input format table. The next-to-last card 
specifies the time factors: the initial time, the final time, 
the integration step size and the frequency of output are 
given in an 8E10.0 format. The last card enumerates the 
variables to be plotted versus time. 
Here some specifics regarding the time specifications 
and the variables to be plotted must be remembered. 
(1) Common sense must be used when selecting the initial 
and final time. Intelligent guesses should be made based on 


experience and the system dynamics. 
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(2) The integration step size is also related to the system 
dynamics. It should be small enough to give a precise solu- 
tion but not excessively small as to increase the computing 
time unnecessarily. As a rule of thumb one can start by 


letting the integration step size be 


(final time) - (initial time) 


= 1000. 


(3) The frequency of output (FREQ) determines both the 
number of points to be plotted in the total time interval and 
the physical dimension of the graph. The formula to determine 
the value of FREQ is 


mme Gime) = (initial time) 


DO = (integration step size) (number of points to be plotted 


= (No. of time steps)/ (No. of points plotted) 


where the number of points to be plotted must always be less than 
or equal to 100. Equivalently one can say that the plotting 


is constrained by the equation 


preo > (final time) - (initial time) 
— (integration step size) (100) 


This relationship is very important. It restricts the user 
but also permits him to establish in advance the number of 
points to be plotted per curve and the scaling of the time 
axis. This is illustrated by the following example. 


Assume that the initial time is 0.0, the final time 


.is 10.0 and the step size is 0.005. What value should be 
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used for the frequency of output, FREQ? Using the rule 


stated above, 


Wl HEN 


OO. 20 


FREQ > 
Thus the frequency of output must be greater than or equal to 
twenty. Expecting a moderately oscillating time response, a 
"number of points to be plotted" equal to fifty is decided 


upon. Thus, 


22070509) 
giving a sampling interval (S.I.) 
S.I. = (FREQ) (Step Size) = (40)(.005) = 0.2 


In summary, for this example, setting FREQ equal to 40 would 
give an output of 50 points, each 0.2 seconds apart between 
the initial value TI = 0.0 and the final value TF = 10.0 
seconds. 

Note that the physical dimension of the graph is 
directly proportional to the number of points to be plotted. 
Fifty points usually gives a good drawing and is suggested as 
starting value. 

(4) Approximate equations for the graph dimension are pre- 
sented as extra information only. These do not help to solve 


the problem but give an idea of what to expect: 
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dependent variable or y axis = 36 cm (fixed) 


independent variable or t axis = (.318)x (number of points) cm 


The last card of the data deck indicates what dependent varia- 
bles are to be plotted. A maximum of eight graphs can be 
output for every program run. If tabular output only are 
desired, the last card is left blank. The variables for which 
time responses are to be drawn are specified by giving the 


symbol that corresponds to the desired variable: 


Symbol Variable to be plotted Symbol Variable to be plotted 
L xl (t) 8 X8 (6) 
2 x2 (t) 9 x9(t) 
3 3:3 VE) S EEGEN 
4 x4 (t) R error signal 
5 2 E) U controller input 
6 x6 (t) x output 
7 aay It? R forcing input 


Table VIII - Symbol Indicating Variables 
to be Plotted by GTRESP 


where the error signal is defined as 


e(t) = r(t) - y(t) 


All the above is summarized by the following table: 


EN 





Entry Input Description Format Columns Used 


1 Problem identification, order 5A4, I2 1-20, 21-22 
of the system (N <10) 

2 Plant matrix A(N xN) (one row 8E10.0 1-10, 11-20, 
per card if N < 8 or one row Seer 
per two cards if N > 8) 

3 Distribution matrix zu (1 x N) 8E10.0 1-10, 11-20, 
(one card if N < 8 or two etc. 


cards if N > 8) ` 


4 Output vector c (1 xN) (one 8E10.0 1-105 11-205 
Garay if N < 8 or two cards etc. 
if N > 8) 

5 feedback coefficients k., Kos 8E10.0 1-10, 11-20, 
k One card if N 8 etc. 
or two cards if N > 8) 

6 Controller gain K 8E10.0 1-10 

7 Initial condition x, (to), 8E10.0 1-10, 11-20, 
X (to), IR.) x(t )” (on One etc. 
card if N « 8 or Lwo cards 
iE N > 8) 

8 Initial time TI, 8E10.0 1-10, 
final time TF, 11-20, 
step size DT, AO 
frequency of output FREQ 31-40. 

9' Any of the following symbols 8Al lee 456,78 


in any of the first eight 
columns of the card (maximum 
of 8): 


eet, lap2,3,4,5,6,/,8,9,A 
Table IX - Input Format Table for GTRESP 


Do Output 
i ee T T 
The problem identification, A, b , c, kK , K, x(t): 
the initial time TI, the final time TF, the integration step 
size DT and the frequency of output FREQ are printed out for 


future reference. Then the tabular output of all the state 
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variables together with the control input u(t) and the output 
y(t) are listed versus time. Finally, the graphical outputs 
are given. As mentioned earlier, one graph is produced for 
each selected variable. At the end of the run, a compact 
solution is presented by plotting all the curves on a single 
graph. 

c. Example 


An uncompensated system is described by 


x] (t) = X» (Œ) 


x, (t) a(t) 

The system is compensated by feeding back both 
states and the graphical time response is to be obtained for 
initial condition only. The initial conditions are x4 (0) = 10.0 
and x, (0) = 0.0. The controller gain equals 1.6. 

The following diagram represents the complete 


system: 


r(t)* 





x, (t) 


Fig 3-16 Feedback System for GTRESP Test 
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Since only the time response to initial conditions is required 


for the problem, r(t) is set equal to zero. The system order 


is two and 


The 


ip 
Il 
fo ol 
¡SA 
¿O 
H 
Il 
o 


TS C = [0 0] 
The feedback coefficient matrix pA = Íl. quen] 
The controller gain K = 1.6 
The initial conditions x (0) = {10. 0] 


From the dynamics of the system, a final time of 10 
seconds is chosen. 


An integration step size of 0.02 is sufficiently small. 


time equations imply that 


(10 - 0) 


(1) EREQ 702770100) 


Wt 
LP 


(2) 50 points are chosen to cover the ten second 
interval so 
_ (10 - 0) E 
NEE S 
(3) One value is going to be plotted every (DT)x(FREQ) 
or 0.2 second. 
(4) The estimated dimensions of the graph can be evaluated 


as 


94 





The variables 


All the above 


dependent variable axis = 36 cm 
independent variable axis = (.318) (50) = 15.9 cm 
to be plotted are u(t), x, (t) and x,(t). 


are entered as specified on the input format 


table and the subroutine RFIND(T,R). The complete computer 


cards set up is then: 


// (standard OS JOB card) ,TIME=2 


// EXEC LINCONF 


//FORT.SYSIN,DD,* 


SUBROUTINE RFIND(T,R) 


R=0.0 
RETURN 
END 
/* 
EDLINK.SYSIN,DD,* 


AAINCLUDE (GTRESP) 


. ENTRY ,GTRESP 


SC? 


USO .SYSIN,DD,* 


GTRESP TEST 


1.0 


02 
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12U 
/* 


Results are shown in Fig. 3-17A-E 


D. MODERN CONTROL SUBPROGRAMS 

The following set of subprograms may be used to analyze 
and design linear feedback control systems which are to achieve 
a specified closed-loop transfer function. 

This group of nine subprograms consists of: the supporting 
subprograms OBSERV, CONOBS, SENSIT which provide the user 
with a means of checking the observability and controllability 
of a system and its sensitivity to parameter variations; the 
subprograms STVAR, LUEN and SERCOM which help design optimal 
linear control systems with complete or incomplete state 
measurements; RICATI and KALMAN which find the feedback and 
control gains necessary to optimize a given function either 
for continuous or discrete systems; finally, MIMO which is a 
computer aided technique to determine feedback control laws 
for multiple-input multiple-output systems where the number 
of inputs equals the number of outputs. 

The subprograms SENSIT, KALMAN and RICATI are normally 
Class B subprograms and require a "TIME = 2" specification on 
the JOB card. All others are Class A. Except for KALMAN 
which must be operated using Mode Two, all the subprograms 


are accessed under Mode One or Mode Three. 
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l. Observability (OBSERV) 


This subprogram determines the observability index of 


the linear, time invariant, Nth order system 


x(t) 


li 
Up 


BIERENG) 


y (t) 


ii 
UC) 


x (t) 


The observability index r of the above system is 
defined [4] as the smallest positive integer for which the 
matrix [C, arc, Ee (A) EEN has rank N. 
IAE 
The problem identification, the order of the system 
Ene number of rows of the C matrix are entered on the first 


Alla deck card. Then the A matrix is presented one row at 


a time followed by the C matrix, also one row at a time. 


Entry Input Description Format Columns Used 
d Problem identification, 5A4 1-20, 21-22, 
order of the system (N < 10), 2 23-24 
number of rows of C (M < 10) 
2 A (N xN) matrix (one row per 8F10.3 1-10, 11-20, 
card if N < 8; one row per etc. 


two cards if N > 8) 


3 C (MxN) matrix (one row 8F10.3 1-10, 11-20, 
per card if M < 8; one Sacer 
row per two cards if M > 8) 


Table X - Input Format Table for OBSERV 
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bi. Output 
The problem identification, and the A and the C 
matrices are listed for reference. Then either "(A,C) is 
unobservable" is printed or the observability index is given. 
(If the observability index equals N the number of states, 
the system is completely observable.) 
c. Example 
The following set of matrices are to be checked 


for observability condition 


-1 -2 -2 
(1) = 0 -1 1l and C = [1 1 0] 
1 0 -] 
2 1 0 


(2) A and 


¿ 
Il 
eo 
N 
RP 
UC) 
i 
(Oz ol 
N p 
e 


Here, both (l) and (2) are solved in the same run, placing the 
data decks one on top of the other giving: 

// (standard OS JOB card) 

// EXEC LINCON 

//LINK.SYSIN,DD,* 

xx INCLUDE SYSLIB(OBSERV) 

p* 

M/GO.SYSIN _DD_* 


OBSERV TEST ONE 0301 
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0.0 24.0 1.0 
1.0 0.0 =l. 0 
1.0 1.0 0.0 


OBSERV TEST TWO 0302 


2.0 1.0 0.0 

0.0 2.0 ISO 

0.0 0.0 Ss? 

0.0 10 SmO 

0.0 20 4.0 ` 
/* 


and the solution is shown in Fig. 3-18. 
2. Controllability and Observability (CONOBS) 
The subprogram is a modified version Of OBSERV. TIE 
is used to determine the observability index and check the 


controllability of a linear, time-invariant control system of 


the form 
x(t) = Ax(t) + B u(t) 
y(t) = C x(t) 
a. Input 


The input is the same as for OBSERV except that 
the B^ matrix must be included. The input deck starts with 
the problem identification card which also contains the system 


order, N, the number of rows of ENDS and the number of outputs, 


104 





ROBLEM IDENTIFICATIONS ~~ OnSeRV TEST CNF 
*9*$9*9*690€99*99-*09]9*****529*62»5$*92*5$60*5*5**9***t€t11zz15 19599 
THE A MATRIX 
-1.90CGCOOE 00 - 
(LGäecoruas og ` 
TEE C MATRIX 
^ 1.00C0COO0E 00 1.000CCOCE 00 0.0 
?*5*$99$9*9294*5*5**49$99929**6*5**5*95595552560**5*x*92**9* 


OBSERVARIL ITY INCEX 3 


E 


)HSERVARILITY TMCEA CALLULATICA 
PRORLEM [OENTIFICATICN CBSERV TEST TWO 
$64969995»99399t«2999599992999922294922* t192*495 


TFE A "ATATX 


2.0000C0JE 00 De 09 n.0 
0.0 2 e CCCCCCCE 00 1.0C00000E€ 00 
0.9 0.0 2.0000000F 00 
THE C MATRIX 
0.9 1.CCCCCCCE 09 3.0000000€ 00 
0.0 2-CCCCCOCE 00 4.0000N00F 00 
dtd dt ttt 


(aC) IS UNCBSEPVABLE 


Figure 3-18 Observability Subprogram Tests 
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M. Next the A matrix (NxN), the SE matrix (LxN) and the 


C matrix (Mx N) are entered one row at a time using an 8F10.4 


format. 
entry Input Description Format 
T Problem identification, 5A4, 


system order (N < 10), 312 
number of rows of B (L < 10), 
number of outputs (M < 10) 


2 A (NxN) matrix (one row per 8F10.4 
card if N < 8; one row per 
two cards if N > 8) 


3 p (LxN) matrix (one row 8F10.4 


poer card TEN < 8; one row 
per two cards if N > 8) 


i C (MxN) matrix (one row 8F10.4 
per card if N < 8; one 
row per two cards if N > 8) 


Table XI - Input Format Table for CONOBS 


b. Output 


Columns Used 


1-20, 21=22; 
23-24, 25-26 


210,120, 
21-30, etc. 


1-10, 11-20, 
21-30, etc. 


1-10, 11-20, 
21-30, etc. 


The problem identification and all three matrices 


are output for reference. Then two sentences are printed 


indicating whether or not the (A,C) system is observable 


and the (A,B) system is controllable. 


c. Example 


The following systems are to be tested for observa- 


bility and controllability: 
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(1) x(t) = 0 -l 0 x(t) + 0 0 uot) 
-3 -4 -2 1 0 
y(t) = x, (t) 
-2 0 Ji 
(2) x = 
x(t) 0 -1 l x(t) T 0 0 u(t) 
-3 0 -2 1 0 
y(t) = x4 (t) 


Here again, both solutions are obtained in a single run using 
one set of control cards before the two consecutive data 
decks. 


Ber (1), 


22 0 1 
NOS = 0 21 0 
-3 -4 -2 
3 Q i 
pal <3) = 
= 1 0 0 
c(1x3) - (1 0 0l 
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Thus system order N 3, number of rows of St EE 
number of outputs M = 1. 


For (2) , 


-2 0 1 
Bere = 0 -1 d 
-3 0 -2 
0 di 
B(2x3) = 
~ l A 0 
C(1 x 3) = (1 0 0] 


and the system is of order N = 3, with L = 2 and M = 1. 
The control cards and data decks are then as follow: 

// (standard OS JOB card) 

/ / ^ EXEC ALINCON 

//LINK.SYSINADD, * 

n INCLUDE SYSLIB(CONOBS) 

/* 

//G0.SYSIN-ADD.* 


CONOBS TEST ONE 030201 


=2.0 0.0 1.0 
0.0 =O 0.0 
ES.0 -4.0 22m 
0.0 0.0 l.0 
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1.0 0.0 0.0 


CONOBS TEST TWO 030201 


=2.0 0.0 1.0 
0.0 =1.0 1.0 
9. 0 0.0 SE 
0.0 0.0 1.0 
1.0 0.0 0.0 
1.0 0.0 0.0 
/* 


The results are presented in Fig. 3-19. 
3. Sensitivity Analysis (SENSIT) 
This subprogram is used to obtain the root locus of 
the closed-loop poles of the (single-input single output) 


linear control system 


x(t) A x(t) + b u(t) 


u(t) G{r(t) - k'x(t)], where G is a scalar, 
as a single element of the plant matrix A,or the control 
vector b, or the feedback coefficients matrix k*, or the con- 
troller gain K varies between some specified values. As 
already mentioned, the subprogram studies the effect of a 
single parameter variation and plots the result. If, for 

the same system, it is desired to consider more than one 


parameter variation the user indicates his choices by providing 
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ORSE«4VAFILITY, CCHT-CLLZEILIT 
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Figure 3-19  Controllability and Observability Subprogram Tests 
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one option card per element to be varied and the computer 
cempletes one root locus for each parameter. The end of a 
problem is indicated by a blank card. After that card, a 
data deck pertaining to other systems may be included if 
desired. 

Execution times for this subprogram are normally more 
than 20 seconds. Thus TIME = 2 should be specified on the JOB 
card. Mode One or Mode Three is to be used to access the 
subroutines. 

a. Input 

The problem identification and the system order 
(N < 10) are presented on the first card. Next the plant 
matrix A (NxN) and the b? eN matrix are entered, followed 
k 


by the feedback coefficients k k and the controller 


1 tx 
gain G. Then the option card is given, indicating the element 
to be varied, the number of parameter values to be used, and 
the minimum and the maximum values of that parameter. This 
card, with the proper modification, is repeated once for each 


element to be varied. Finally a blank card indicates the end 


of the problem. 


LII 





Entry 


5 


6 
(repeat 
this 


Input Description 


Problem identification, 
order of the system (N S B) 


A (NxN) matrix (one row per 
card if N < 8; one row per 
two cards if N > 8) 

pi (lxN) matrix (one card 

if N < 8; on two cards if 

N > 8) 

T Bea e 
k (1 xN) coefficients matrix 
(one card if N < 8; two cards 
if N > 8) S 


Controller gain G 


(1) element to be varied 
(letter A if the 


entry once element is part of A, 


for each letter B if the 
parameter element is part of b, 
which is K if the element is. 
to be one of the feedback 
varied) coefficients, G if the 


element is the 
controller gain) 


(2) row number of the element 
duerrom A, b, or k. 
Otherwise set equal to 
00. 


(3) column number of the 
element if from A. 
Otherwise set equal 
to 00. 


(4) number of parameter values 
to be used. 


(5) minimum value of the 
parameter. 


(6) maximum value of the 
parameter. 


blank card (this indicates 
the end of the problem) 


Table XII - Input Format Table for SENSIT 
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Format 


5A4 
12 


8F10.3 


SE TOSS 


El 


82103 


Al 


I2, 


I5, 


210.37 


Faro 3 


(blank) 


Columns Used 


l-20, 21-22 


l-10, 11-20, 
etc. 


esM 2E 
etc. 


al 207 
Ee 


is O 


6-10 


a2 0 


ZU 


(blank) 





The user must be very careful while preparing the 
option cards. The following example can best illustrate the 
procedure. Suppose it is desired to get the root locus of 
the poles of a closed-loop system as the parameter a34 varies 
from 0.0 to 100.0. A "number of parameter values to be used" 


of 20 is selected giving the following option card: 


colum 123456789 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 
data ARORO 0 0 . 0 1 0 0 


If for the same problem, it is also desired to study the 
variation of the closed-loop poles as b, varies from 0.0 to 
100.0 with a "number of parameter values to be used" of 10 and 
also as G varies from -1600. to -1200. with a "number of 


parameter values to be used" of 25, then the two added option 


cards would be: 


tu 1234567891011 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 


data B0300000100. 0 100. 
data G000000025 - 1 6 00. - 1200. 
b. Output 


The problem identification, the A, po, k? matrices 
and the gain value G are listed first. Then the first element 
to be varied and its minimum and maximum values are printed, 
followed by each parameter value and the closed-loop poles 
associated with it. Finally the root locus plot is given for 


each element to be varied. 


ITS 





Note that the "number of parameter values to be 
used" should be kept small. Since the values of the roots 
are calculated and printed for each parameter value, 100 
values should be regarded as a practical maximum. 

C. Example 

The stability of the following system is to be 

investigated under gain variation and the effect of the non- 


perfect integrator => looked at. 


T(t)+ 





Fig 3-20 Control System for SENSIT Test 


First the state equations 


EQ - Sx (t) + x (t) 
x, (t) = -2x, (t) + X, (t) 
E) (E) + u(t) 

u(t) = Glr(t) - x, (t)] 


are written, giving the following data: 
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Order of the system: N = 3 


-€ l 0 
Ewe 0 -2 1 
0 =; 
T 
b° = [0 0 1] 
eege 0 0] 


The gain G is to be varied from 0 to 10 with e set to zero, 
and a "number of values to be used" of 20 seems reasonable. 
The value of e is also to be varied with G set to its nominal 
value of 1.0. First e is set equal to zero in the A matrix. 
A range of variation of 0 to l and a "number of values to be 
used" of 25 are selected. 

The control and data cards are: 

// (standard OS JOB card), TIME-2 

// »EXECALINCON 

WABENK.SYSINADD _* 


^^INCLUDE^SYSLIB (SENSIT) 


/* 

E60. SYSIN.ADD.A* 

SENSIT TEST 03 
0.0 1.0 0.0 
0.0 22:0 1.0 
0.0 0.0 D) 
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G0000000200.0 10.0 
A010100025 0.0 1.0 
(blank card) 
/* 
The results are presented in Fig. 3-20(A-D). 

4. State Variable Feedback (STVAR) 

This subprogram is a very powerful aid for design and 

analysis of any single-input single-output linear, time- 


invariant system represented by the states equations 


x(t) = A x(t) + b u(t) 
T 

BE EE) - k x(t)i 

v(t) = ¢ x(t) 


It permits one to find internal transfer functions 
of the form X: (s) /U(s), the plant transfer function Y(s)/U(s), 
the closed-loop transfer function Y(s)/R(s) and the equivalent 
feedback transfer function Hog (S). In addition, this subpro- 
gram calculates the controller gain and the feedback coeffi- 
cients necessary to achieve a specified closed-loop transfer 
function. It is to be run as a Class A job and is accessible 


under Mode One or Mode Three. 
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The procedure is not very complex, but 
requires understanding. All the information needed by 
the user to solve state variable feedback problems is 
presented in the following paragraphs. However, the tNesryson 
which the subprogram is based is not given. The user who 
wishes to learn more about it should refer to the texts by 
Schultz and Melsa [5], Melsa and Jones [1], Eveleigh [6] or 
others. 

a. nbi 

(1) Basic Cards 
As usual the problem identification and the 

system order are given on the first data card, followed by 
the plant matrix A (Nx N) and the transpose control vectors 
b? (1xN). From this input (which is always reauired) the 
subprogram verifies the controllability of the system. Three 
possible controllability conditions may be found by the 
computer. One, the system is completely controllable and no 
special message is printed. Two, the system is numerically 
uncontrollable. In other words, it is theoretically controlla- 


ble but uncontrollable in a numerical sense. This situation 


arises when the controllability matrix 


E = [B AB AB MM X5 


ke, gd Pe Ma mo MI v2 


cannot be accurately inverted using the programmed algorithm. 
The matrix and its calculated inverse are then multiplied 


together and checked against the identity matrix to provide a 
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measure of the uncontrollability of the plant. If the des- 
cribed condition occurs, the message "plant is numerically 
uncontrollable" is given accompanied by "MAX. DEVIATION=number", 
where "number" is the value of the deviation from the identity 
matrix. Reference l states that a maximum deviation value 


> has been found to indicate diffi- 


larger than on to 10° 
culty. The last controllability condition is "the system is 
uncontrollable", and is indicated as such. Note that even if 
the plant is determined to be uncontrollable, the computer 
solves the problem and presents the results. The option of 
accepting or rejecting the solution is the designer's preroga- 
tive. 
(2) Open-Loop Cards 

The next input cards specify which open-loop 
transfer functions are to be computed. These cards need not 
be provided if no internal transfer function is desired. The 
way to identify the internal transfer functions to be output 
is by using ficticious c matrices. The following example 
demonstrates the procedure. Suppose the internal transfer 
functions X, (s) /U(s) and X1 (s) /X, (s) are desired for à fourth 
order system. Since only Xi (s) /U(s) type of transfer functions 
are computed by the subprogram, Xj (s)/U(s), X) (s) /U(s) and 
X, (s) /U(s) are requested and the user then only needs to 
divide X1 (s) /U(s) by X,(s)/U(s) to obtain X1 (s)/X,(s). The 


fictitious c, matrices to be provided as input are then: 


f 
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X, (s) 


(1) for U(s) ' 
ce = Il 0 0 0] 

X, (s) 

(2) for UNS) 
Ce = [0 1l 0 0] 

X,(s) 

( 3) for U(s) ' 
Ge = [0 0 0 1] 


Following these cards, the real output matrix c and a null 
matrix DEIN) must be entered. The real c matrix is used 
to compute Geif and correctly solve the rest of the problem. 
The 0 matrix is necessary to indicate the end of open-loop 
calculations. 

(3) Closed-Loop Cards 

Finally the closed-loop input data are given. 
Here again the user has a choice among three options. 

The first of these closed-loop computations 
is for analysis only. This choice is indicated by an option 
card on which the letter A is printed in column one. Following 
this card, the feedforward gain K and the feedback coefficient 
matrix E are given as specified on the input format table. 
From this input, the subprogram determines the closed-loop 
characteristic polynomial and the numerator of the equivalent 


feedback transfer function (both the factored and unfactored 
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forms). From these, the block diagram shown in Fig. 3-21 
can be drawn where u = Y(s)/U(s) and Hog (8) is the 


equivalent feedback transfer function. 





Fig 3-21 ea Form Block Diagram 


The other two closed- loop computations are 

for design purposes. They are used to calculate 

the controller gain and the feedback coefficients necessary 
to achieve a given closed-loop Bo polynomial. 
This desired polynomial is the denominator of Y(s)/R(s) and 
must agree with the system order. If the characteristic 
polynomial is to be entered in P form, an option card with 
the letter P in column one is presented followed by one (if 
E a: ker n 2.8) cards containing the coefficients in 
ascending order. The coefficient of the highest degree term 
must always be 1.0 and may be entered as ten blank spaces. 
On the other hand, if it is more convenient to present it in 
factored form, the option card has the letter F in the first 
column and the next cards give the real and imaginary parts 


of the root using a 2E10.0 format. 
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Since a user may very well wish to obtain 
the closed-loop computations for many different characteristic 
polynomials or try out several values of feedback or feed- 
forward gains, the subprogram allows one to ask for as many 
closed-loop computations as desired by placing the input 
cards one on top of the other. 

(4) Problem Termination Card 

The last card must be blank. It indicates 
the end of the problem and must always be present, whether 
or not the closed-loop portion 35 included. The following 


format table conveniently summarizes all the above. 


Entry Input Description Format Columns Used 

I Problem identification, 5A4, I2 1-20, 21-22 
(Basic) system order (N < 10) 

2 Plant matrix A (NxN) (one 8E10.0 1-10, 11-20, 
(Basic) row per card for N < 8; one etc. 

row per two cards for N > 8) 

3 Control vector b* (1xN) TOO 

(Basic) (one card for N < 8; two etc. 


cards for N > 8) 


4 Cf Guo (one card for 8E10.0 1-10, 11-20, 
epen- "N «79; two cards for N > 8) etc. 
loop) (repeat if several fictitious 

matrices) 

5 Output matrix c (1xN) (one 8E10.0 1-10, 11-20, 
(open- card for N < 8; two cards etc. 
loop) för > 8) 

6 Null matrix 0 (1xN) (one 8E10.0 all 
(end of blank card for N < 8; two 
Open blank cards for N > 8) 
loop 
part) 
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Entry 


7 
Analy- 
sis 


8 
Analy- 
sis 


9 
Analy- 
sis 


10 
Design 
option; 
unfactored 

form 


EL 
Design 
option; 
unfac- 
tored 

form 


12 
Design 
option; 
factored 
form 


13 
Design 
option; 
factored 

form 


14 


Input Description Format Columns Used 
Letter A in column one Al 1 
Feedforward gain 8E10.0 l-10 
Feedback coefficient matrix 8E10.0 MSN 207 
kT (1xN) (on one card for etc. 

N < 8; two cards for N > 8) 

Letter P in column 1 Al 1 

Desired characteristic poly- 8E10.0 1=10, 11-20, 
nomial coefficients (one one etc. 

cod N < 8; two cards if 

N > 8). See p. (31) for 

details 

Letter F in column 1 Al 1l 

Desired characteristic poly- 8E10.0 IO TI-20. 
nomial roots (one per card, 

real part followed by 

imaginary part. See p. (32) 

for details 

Blank card (indicates the 8E10.0 blank 


end of this problem) 
Table XIII - Input Format Table for STVAR 


Note that entry (4) must be included if no internal transfer 


function is desired. 


The same also applies to entries (7-8-9) 


if analysis option is not desired and (10-11) and/or (12-13) 


if no design option is taken. 


ES 


Output 


The problem identification and the A and DE 


matrices are given and, if applicable, a numerically or 
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completely uncontrollable situation is indicated. Next the 
open-loop calculations are presented. The denominator 
coefficients in ascendina powers of s and the roots of the 
denominator polynomial are listed at the beginning of the 
section. Then, if requested, each ficticious Ce matrix 
followed by the numerator of the corresponding transfer func- 
tion is printed. The last output of this section is the C 
matrix and the numerator of the plant transfer function. 
The user is reminded that the ficticious Ce matrices indi- 
cate which X; (s) /U(s) is computed while the c matrix specifies 
the real output y(t) which is used to calculate Y(s)/U(s). 

The last section of the printout concerns the 
closed-loop calculations. If the analysis mode was selected, 
"KEY=A" is printed followed by the numerator of the equiva- 
lent feedback transfer function, DE both in factored and 
unfactored forms. Note that the complete Sl is obtained 
by taking the "numerator of Hog (given in the closed-loop 
calculations) and dividing it by the numerator associated with 
the real c matrix (given as the last part of the open-loop 
calculations). Next the feedback coefficients and the gain 
are listed for reference and the computed closed-loop char- 
acteristic polynomial and its roots are given. 

If computations of the feedback coefficients and 
the feedforward gain to achieve a desired closed-loop char- 


acteristic polynomial was requested, the computer output 


shows "KEY=P" or "KEY=F", depending on the design mode 
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selected. Then, as for the analysis mode, the numerator of 
the Hog (5) is given, followed by the feedback coefficient 
matrix B and the feedforward gain K. Here it must be pointed 
out that the subprogram calculates the gain K so zero steady- 
State error results froma step input. A designer who wishes 
other conditions may rescale K and k? appropriately by hand. 
For example, suppose it is desired to have the controller 

gain K = Ky but the computer output shows that K = Ko with 

l' K, and K,. The procedure is 
then to modify the results by setting 


the feedback coefficients k 


and setting 


K = —[k K Kew 


This does not change Y(s)/R(s) and satisfies the condition 
K = K}. Finally a parameter called "maximum normalized 
error" is associated with each closed-loop calculation. The 
value of this parameter indicates the exactitude with which 
the problem was solved by the computer. This number can help 
to determine the validity of a solution, especially when 
numerical uncontrollability was encountered to start with. 

c. Example 


Eveleigh [6] presents the ideas of design of 


control systems using state-variable feedback and works out 
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two examples, the first of which is solved here by the 
computer method described previously. The problem can be 


stated as follows: given the plant transfer function 


10 


3p) 7 s(s*0(s*3)' 


find each state feedback gain and the feedforward path gain 


necessary to achieve the closed loop transfer function, 


BHL m > 10 


S +45? +9s +10 


The first step of the procedure is to get the 


state variable representation of the system. The following 


signal flow graph may be obtained: 





Fig 3-22 Signal Flow Graph for STVAR Test 


By inspection, 


0 l 0 
A = 0 zu 2 
0 0 = 
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Normally the user first runs the subprogram for 
open-loop calculations. Then he either uses it for analysis 
or for design. To illustrate all the possibilities, the 
subprogram was applied to solve the same problem using all 
of the different modes. 

(1) Open-Loop Test 

For the case at hand, assume the solution is 
to include the internal transfer function X, (s) /U(s). Thus 


the input data requires a ficticious c, matrix to be added, 


f 


l.e. 


The computer card deck Cor this simple 
open-loop test is: 
// (standard OS JOB card) 
// EXEC LINCON 
ELINK .SYSIN.DD,* 
n n INCLUDE,SYSLIB(STVAR) 
/* 
//GO.SYSIN,DD.* 


STVAR OPEN LOOP TEST 03 


130 





0.0 =1.0 2.0 
0.0 0.0 E 
0.0 0.0 5.0 
0.0 1.0 0.0 
1.0 0.0 0.0 


(blank card) 
(blank card) 
par 

The first blank card is a null matrix 0 (1x3) 
that indicates the end of open-loop calculations while the 
second blank card indicates the end of the problem. From 


the results shown in Fig. 3-23, 





SEN, > 10s | 10 

Us) | s(s+1)(s+3) (s+1)(st+3) 
and 

NUS . 10 

Uls) | s(s+1l)(s+3) 


(2) Analysis Test 
To illustrate the analysis computations, the 
feedforward gain K = 1l and the feedback coefficients Ky = 1, 
k- = 0.6 and k_ = 0 were assumed. Again the computer card 


2 3 
deck is given below. 


SE 











STATE VARIABLE FEECBACK 
PROILEM IDENTIFICATION - STVAR OPEN-LNOP TEST 
SSSSHSsSeseStKe AAA ARE 


THF A MATRIX 


0.0 1-0006COOE 00 0. 0 

0.0 -1.0000000€ 00 2.0000000E 00 

Gen 0.0 -3.0000000E 00 
THE B MATR IX 

0.0 0.0 i 5.9000000E 00 


TTT TH CTT OCT TTC CTT ET ECH 
OPEN=LOOP CALCULATICAS 
DENOMINATOR COEFFICIENTS - IN ASCENDING POWERS OF S 


0.0 3.0000U00€ 00 . 4.0000000E 00 1.0000000E 00 
THE ROOTS ARE REAL PART IMAGINARY PART 
-1.CCQ0000E 00 “Oe 
-3.0000000E 00 0.0 
C.C 0.0 
THE C MATRIX T 
0.0 1.0000000€ 00 . 0.0 
" NUMERATOR CCEFFICIENTS - I^ ASCENDING PCWERS OF S 
0.0 1.0000000E 01 
THE AONTS ARE REAL PART IMAGINARY PART 
CO 0.0 
THE C MATRIX sesso 
1.0000000E 00 0.0 0.0 


NUMERATOR COEFFICIENTS — IN ASCENDING POWERS OF S 
1.0000000E Ol 


Figure 3-23 State Variable Feedback - Open-Loop Test 
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// (standard OS JOB card) 
// EXEC, LINCON 
//LINK.SYSIN DD,* 

Wi INCLUDE SYSLIB(STVAR) 
pet 

MAGO. .SYSIN DD * 


STVAR ANALYSIS TEST 03 


0.0 160 0.0 
0.0 Fe) 2.0 
0.0 0.0 25.0 
DO 0.0 3.0 
Be 0.0 0.0 


A 

2.0 

1.0 0.6 0.0 
(blank card) 

p" 


Interpretation of the output reproduced in 


Fig 3-24 gives 


_ 10 + Gs 
eq 10 
and shows that the closed-loop poles are at ~2 and -l+tj2. 
(3) Closed-Loop Test 
Here the subprogram is used for design. 


Suppose that the feedforward gain and the feedback coefficient 
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STATE VARIABLE FEEDBACK 
PROBLEM IDENTIFICATICN - STVAR ANALYSTS TEST 


LOIRA IAEA INCA CAICARA ICAIC 


ThE A MATRIX , 


0.0 1.6C000000€ 00 0.0 

9.0 -1.0C0cCOC* OQ 2.0000000€ 00 
0.0 0.0 ^ -3.0000000E 00 
THE B MATRIX 

0.0 0.0 5.0000000E 009 


ICT CTAA TTT UCI KIT ZC? 
OPEN-LOCP CALCULATIONS 
OEMNIMINATOR COEFFICIENTS - IN ASCENDING POWERS OF S 


0.0 3.,0C0C0CCE OO 4.0000000€ 00 1.C000000€ 00 
THE ROOTS ARE REAL PART IMAGINARY PART 
-1.C000000F 00 0.0 
0.0 0.0 
THE C MATRIX *..... 
1.0000000E 00 0.0 0.0 


NUMERATOR COEFFICIENTS - IN ASCENDING POWERS OF S 
1.0000000E 01 
$959$99$5*59$95925»9599»459595959259$559*9*55995*$*5959595429*5*$5*959 
CLOSEO-LOOP CALCULATIONS 
KEY = A eso 
THE NUMERATCP JF F-ECUIVALENT - I'l ASCEHCING POWERS OF S 


1.0000000E Ol 5.9599590€ 00 0.0 . 
THE ROOTS ARE REAL PAPT IMAGINARY PART 
-1.6666660E 00 0.0 


THe FEEDIACK COEFFICIENTS 
1.0000000€ 00 5.9999556E-01 0.0 
ThE GAIN = 1.0060000€ 029 
THE CLOSED-LOJP CHARACTERISTIC PCLYNOMTAL — IN ASCENDING POWERS OF S 


1 .0000000E OL .99595681E 00 4. 0000000E 00 l.e 0000000E 00 
The ROOTS ARE REAL PART IMAGINARY: PART 
-1.CC00000E 00 -2.0000002E 00 
-1420000000-* 00 2.0000000€ 00 
-2.0000000E 00 0.0 


MAXEMIIM NIKMALIZED ERROR = 1.06F-07 


Figure 3-24 State Variable Feedback - Analysis Test 
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values are to be obtained so the closed-loop characteristic 


3 2 


polynomial is s +4s +9s+10 or, equivalently, the closed- 


loop poles are located at -2 and -1 * j2. 


For illustration, 


calculations are requested for both the P and the F forms 


(even though they are exactly the same). 


and data deck are then: 
// (standard OS JOB card) 
V7 RESEC.LINCOM 

FEINE. SYSIN DD $ 

ER INCLUDE SYSLIB(STVAR) 
/* 

//GO.SYSIN,DD,* 


CLOSED LOOP TEST 03 


0.0 O 0.0 
0.0 zm SÉ 
0.0 0.0 =3 0 
0.0 0.0 DO 
1.0 0.0 0.0 


(blank card) 


E 

n.0 20 

2.0 

P 

10.0 920 4.0 dee 


(blank card) 


A 
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The control cards 





As expected, the results shown in Fig 3-25 
specify a gain of one and feedback coefficients values of 
ki = 1.0, K, = 0.6 and K3 = 0.0. 

Note that all the above calculations could 
have been executed as a single run using the following card 
deck: 

// (standard OS JOB card) 
// EXEC, LINCON 
M/LINK.SYSIN,DD,* 

a a INCLUDE SYSLIB (STVAR) 
/* 


7760 .SYSIN,DD,* 


STVAR TEST 03 
0.0 MEO 0.0 
0.0 zu ZIO 
0.0 0.0 230 
0.0 0.0 S 
0.0 1=0 0.0 
T. O0 0.0 0.0 


(blank card) 


P 


E. 0 
10 0.6 0.0 
(blank card) 
/* 
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STATE vARIADLE FEEFRDACK 
PROBLEM IOENTIFICATICN - CLOSFO-LOOP TEST 


BER SCHSSENTACSATSEE SLR ETCH TKS ESRRESABREKLAKEREKE 


THE A MATRIX 


0.0 1.0C00000% 00 0.0 
0.0 -1.0000000É 00 2.0000000* 00 
0.) Ye 0 -2.0003000€ 00 
THE B MATRIX 
0.0 0.0 5.0000099E 00 
SSS SHASTSSETCSSTEKCSHSSCHSTCA THs Sc sa seas ese sss eee aeeee 
OPEN=LNCP CALCULATIOAS \ 
DENOMINATOR COEFFICIENTS = IN ASCENDING POWERS OF S 
0.0 3.00000U0E 00 .4.0000000€ 00 1.0000000€ 00 
THE ROOTS ARE REAL PART IMAGINARY PART 
-1.0CJ000UE 00 0.0 
-3.0000000É 00 0.0 
0.0 0.0 
THE C MATRIX $*s992 
1.0000C00E 00 0.0 |. 0,0 


NUMESATCR CCEFFICIENTS - IN ASCENDING POWERS OF S 
1.0090000E OL 
4*$$€$*992545*€xs$2$*«9$t*shty»*29g9»£&f?*92&k&t929993923*»*»x*9$*9$95»x*X4$9*99 
CLOSEU-L02P CALCULAT IONS 
KEY = F ease 
THE N'UMERATOP OF F-EQUIVALENT - IN ASCENOING POWERS OF S 


1.000U0Us+E Ol 6.0C00019E NO 0.0 
THE ROCTS ARE RELL PART IMAGINARY PART 


-1.€666660£ 00 0.0 

THe FEECAACK COEFFICIENTS 
1.00u0000€ 00 6.0CGCCC2E-01 0.0 

Tre GAIN a 92-9999 ISBE-OL 

THE CLOSED-LOOP CHARACTERISTIC PCLYNOMIAL - IN ASCENDING POWERS OF S 


$.9939962€ 06 8.9999962E 09 4.00J0000E 00 1. 0000000€ 00 
THE ROOTS ARE REAL PART IMAGINARY PART 

-1.0C00009E 00 -1,9999990F 00 

-1.0009090E 0N 1.9999990E 00 

~2e0COf000S 00 e 


MAXIMUM NORMALTZED ERRCR = 3.19E-07 


KEY = P “>. 
THE NUMERATONR OF HP-EQUIVALENT - IN ASCENDING POWERS OF S 
1.000200*4E Ol 6.0006015< 00 0.0 


THE RCOTS AFE PAST DET PART 


-16656600E 00 
THE FEEDBACK COEFFICIENTS 
1.00C0000€ 00 6.000u002E-01 0.0 
TRE GAIN 2 9.,9999958E- C1 E o A 
THE CLOSED-LOOP CHARACTERISTIC PPLYMOMIAL =. IN ASCENDING POWERS OF S 


9.9999962E 00 8.99995€2E 00 4. 0Cc00000F 00 1. 0000000E 00 
THE ROOTS ARE REAL PAF? IMAGINARY PART 

-1.CCO00N0E 00 -1.9999990€ 00 

-1.0600000* 00 1.9999990E 00 

-2 ..0000000< 00 © 


MAATMUM NORMALIZED ERROR = 3.18E-07 


Figure 3-25 State Variable Feedback - Closed-Loop Test 
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(4) Step Procedure 

The procedure demonstrated through this 
simple example applies for all problems. The steps to be 
taken can be summarized as follows: 

(a) Obtain the state variable representa- 
tion of the system. 

(b) get A, b", and c. 

(c) If necessary, define ficticious Cf 
matrices to compute "internal" transfer functions. 

(d) For analysis, select the feedforward 


k k 


1! Kar rer Ku. 


(e) For design, select the desired closed- 


gain K and the feedback coefficients k 


loop characteristic polynomial or poles to be achieved. 
>. Luenberger Observers (LUEN) | 

The subprogram LUEN is used to design a combined 
observer-controller to achieve a desired closed-loop transfer 
function when some of the states are not accessible. The 
following paragraphs present a detailed description of the 
computer aided design procedure. However the theory of 
Luenberger Observers in the design of linear, time-invariant 
feedback control is not included in the discussion. Users 
who are not familiar with the subject should consult refer- 
ences 4 and 7, or any other relevant textbook before working 
with this subprogram. 

The solution plan is to start from the state variable 
representation Of a linear time-invariant system and recon- 


struct the missing states using an observer. Then, using 
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both measured and estimated states, assign the feedback 
coefficients and gains required to properly control the 
system. The block diagram presented in Fig 3-26 best shows 
what is meant. The plant represented by the state variable 


equations 


x(t) A x(t) + b u(t) 


y (t) 


Il 
O 


x(t) 


must be controllable and observable. Notice that the C 
matrix indicates which state variables are measured. For 
example, a fourth order system with only the states Xo 
and X4 being accessible would yield 


O 
N 
o o 
O H 
e O 
3 Or 


The real output to be controlled, denoted by SUN may 
elther be one of the state variables or a linear combination 
of several of them. The user is to define a desired closed- 
loop transfer function and find what feedback gains would 
normally have to be used to obtain it, assuming all states 
were available. This is done using the subprogram STVAR 
as explained later in the design procedure. 

The subprogram LUEN is then used to calculate all 


the elements necessary to construct the observer and the 
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Bit 


y_(t) 
C 
| cwmenzm | Er | am es 


y (t) 


Fig 3-26 Luenberger Observer Block Diagram 
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controller. The designer only has to specify, in an arbi- 
trary manner, the observer eigenvalues and the necessary 
feedback coefficients previously found by the use of STVAR. 
The computer solution gives all the matrices and gains 
required. Brought together these form the following compen- 


sated system: 


II E) 
E ASS 
u(t) = Klr(t) - g y(t) - h^x(t)] 
y(t) = € x(t) 
where 

x(t) = state vector 

UNE) = input to the plant 

y(t) = output vector 

Et) = system forcing input 

x(t) = estimated state vector 

A = plant matrix (NxN) 

B = distribution matrix (N xl) 
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F = observer eigenvalue matrix 


G AP = Observer gain matrices 

K controller’ gain 

jt E: : 

g = output feedback coefficient matrix 

T a ; 
h = observer feedback coefficient matrix 
€ = output matrix. 


All these elements except for K, which comes from STVAR 
results, are given as output of the subprogram LUEN. The 
four equations defining the compensated system can be easily 
rearranged, as demonstrated in the example which follows, 
to simulate the system by the use of the subprogram GTRESP. 
a. Design Procedure 

The step-by-step design procedure presented here 
contains the essential information to use the program. It 
also summarizes the Luenberger Observers design concepts. 

Seen. 1 

The closed-loop transfer function Y.(s)/R(s) to be achieved 
must be selected and, assuming all states to be measured, we 
solve for the controller gain K and the feedback coefficients 
K 


k k_* This is done using the state variable 


kt DE s est 
feedback subprogram STVAR, which also checks for system 
controllability. It must be kept in mind that the c matrix 


for STVAR is the matrix associated with the real output yo lt) - 


142 





Step 2 


If an acceptable solution resulted from STVAR, the observa- 
bility index must next be determined. This can be done by 
the use of the subprogram OBSERV, or by hand, using 


G = ter atc? (qaTyj?gT ... (ar) ro lor, 


where the observability index r is the minimum integer such 
that the matrix G has rank r. If (A, C) is found to be 
observable, an observer whose order is equal to or greater 
than (r-1) can be designed. 

Beep 3 
The eigenvalues of the observer are selected arbitrarily. 
However, to ensure a unique solution will exist, TEVIS 
necessary to let the eigenvalues of F be different from those 
Of A. The eigenvalues of A were previously calculated by 
STVAR so it should be very easy to choose some appropriate 
roots for the observer. 

Step 4 
Being the input format for LUEN, the data are entered and the 
subprogram executed. The following system is the final 


result: 


x(t) - 


T 


x(t) - b NE) 


y (t) - 


de 


x(t) 
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X (t) - F x(t) + Gi Ware) n G3 u(t) 


u(t) 


Kee = gy (te= h k(t) 


Bep 5 


If desired, the above equations are rearranged using simple, 


although sometimes laborious, matrix manipulation as 


x(t) x(t 
^ = A . + br(t) 
x(t) 8 x(t) S 
MB) = c x(t 

RS 


Note that the above augmented system order is equal to the 
order of the plant, N, plus the order of the observer. The 
complete system is finally simulated by the use of GTRESP 
letting k? equal zero and K ecual to unity. 
DEN SUE 

As usual the data deck begins with the problem 
identification card on which the order of the plant, N, the 
number of measurements M and the order of the observer, (r-1) 
or greater, also appears. Next, the plant matrix A (Nx N), 
the distribution vector b (1xN) and the measured states 
matrix e (Mx N) are given one row at a time. The feedback 


coefficient matrix kl is then entered exactly as given by 
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the state variable feedback subprogram (STVAR) output. 
Finally, the observer eigenvalues, which are different from 
those of the plant, are supplied either in the form of a 
characteristic polynomial, option P, or as the roots of that 
polynomial, option F. The option is specified in column one 
of the first card by writing the letter P or the letter F. 
If option P is selected, the characteristic polynomial coeffi- 
cients are given in the usual ascending order fashion, with 
the highest order coefficient always set equal to 1.0. For 
example, if the characteristic polynomial of a third order 
observer is chosen to be oe Bee the last two data 
deck cards would then be: 

P 

16.0 4.0 5.0 1.0 
On the other hand, if the roots are to be entered as such, 
the letter F is written on the option card followed by the 
observer eigenvalues presented in the usual manner. For 
example, if the observer poles are -2, -2, -1+3, -1-3 the 


cards would then be 


1.0 1.0 
Note the sign inversion and the fact that only the complex 


root with the positive imaginary part is entered. 


The following input format table summarizes 


the above. 
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ENCTYy 


Input Description 


Problem identification, 
order of the plant N < 10, 
dimension of the output 
vector M, order of the 
observer L  (r-1) 


Plant matrix A (NxN) (one 
row per card for N < 8; 

one row per two cards for 

N > 8) 

Distribution matrix pT (1 x N) 
(one card if N < 8; two 
cards if N > 8) 


Measurement matrix C (MxN) 
(one row per card for 

N < 8; one row per two cards 
for N > 8) 


Feedback coefficient matrix 
k” (1xN) (on one card if 
N < 8; two cards if N > 8) 


Letter F (if observer eigen- 
values are to be entered as 
Beats) or letter P (if 
observer eigenvalues are to 
be supplied by giving a 
characteristic polynomial) 


Entered the observer eigen- 
values as specified on the 
Previous card. (If option F, 
enter the roots real and 
imaginary parts; if option P, 
give the characteristic 
polynomial coefficients 

in ascending order). 


Table XIV - Input Format Table for LUEN 


REGELE 


Format 


5A4, 
3I2 


851073 


SED 


GETO.S 


SIENTO 3 


Al 


8F10.3 


Columns Used 


zum 21-22, 
23223, 223-206 


1-10, 11-20, 
21-30, etc. 


1-10, 11-20, 
21-30, etc. 


LL 2 0 
21-2305 6tc. 


leo, 4-20, 
21-30, etc. 


1-10, 11-20, 
21-30, etc. 


bn T 
The problem identification followed by the A, b 


and C matrices, the desired feedback coefficients and the 
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observer eigenvalues, both in factored and unfactored form, 
are presented for reference. The observer and controller 
elements are printed next as the F matrix, the G, matrix, 


the G matrix, the output feedback coefficients gi and the 


2 
compensator feedback coefficients hl. 

The complete solution of a problem should also 
include the results from the subprograms STVAR, OBSERV and, 
if a simulation is oerformed, GTRESP. 

d. Example 

The example presented by Eveleigh [6] pp. 357-360 

was slightly rearranged and the state variables X, and Xy 


were assumed to be inaccessible. The signal flow graph for 


the uncompensated system is then: 





Fig 3-27 Uncompensated System for LUEN Test 


where x is the controlled output and x4 (t) and x» (t) 


are the measured states. 
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From the diagram, the system matrices are: 





1 0 0 
0 1 0 
A = 
` 0 0 i 
-15  -23  -$| 
T 
m - [D 0 0 1] 
c = E gar 


The solution presented next utilizes the design procedure of 





part a. 
Seep 1l 
The closed-loop transfer function to be achieved is chosen 
to be: 
BER _ I 
RE 


| l 
SES 2) (SL +32) (s + 1-32) 


The controller gain K and the feedback coefficients required 
are found by the use of the subprogram STVAR for which the 


control cards and data deck are: 
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STATE VARIABLE FEFCBACK 
PROBLEM IOENTIFICATION - STVAR FOR LUEN TEST 


$$99***2€6$4$59999$96:,$994$95959*8:5:€99»5929999$99$529522»259255929609 


The A MATRIX 


0.0 1.0000000E 00 0.0 0.0 
9:0 0:0 d OM ce 1290000008 00 
0.0 -1.5000000E OL -2.3000000€ 01 -9.0000000€ 00 
THE 8 MATRIX 
0.0 0.0 0.0 1.0000000€ 00 


$959$*229999949€99999$969299»999269 $22 9095949 $5999399 
QPEN-LOOP CALCULATIONS 
DENOMINATOR COEFFICIENTS - IN ASCENOING POWERS OF S 


» 


0.0 1e-5COCOCOE Ol 2«3C00000F O1 9.0000000E 00 
THE ROOTS ARE REAL PART MA Y l 
=3.0C00000€ 00 De RART 
-9.9999994E-01l 0.0 
0.0 0.0 
THE C MATRIX 10089 
2.0000000€ Ol 1.0000000E 01 0.0 0.0 


NUMERATOR COEFFICIENTS = IN ASCENDING POWERS OF S 
2.0000000€ 0l 1.0000000E Ol 


THE ROOTS ZRE REAL PART 
-2.0C00000= 00 


seta std 0 0 TLT KR DLR 
CLOSED-109P CALCULATICNS 
KEY = F *s*9?9* 
THE NUMERATOR OF H-EQUIVALENT - IN ASCENOING POWERS OF S 
2.0000000E 01 1.3COCCCOC Ol -6.0000000€ 00 


IMAGINARY PART 
0.0 


-3.0000000€ 00 


ThE ROOTS ARE REAL PART IMAGINARY PART 
-1.0396233= 00 0.0 
3.2062893€ 00 0.0 


THE FEEDBACK COEFFICIENTS : 
2.0000000E Ol 1.3000000E Ol 
THE GAIN = 1.0000U900E 00 

THE CLOSED-LOOP CHARACTERISTIC POLYNOMIAL - IN ASCFNOING POWERS OF S 
2.00C0000E Ol 2eS00CCCOE Ol 1.7000000E 01 6.0000000E 00 


-6.0000000E 00 


THE ROOTS ARE REAL PART IMAGINARY PART 
-1.0009000E 00 -2.0000000E 00 
-1.0000002t* 00 2.0000000€ 00 
-2.0004845F 00 0.0 
-1.9995155E 00 0.0 

MAXIMUM NJRMALIZED ERROR = 0.0 


Figure 3-28 STVAR Results for LUEN Test 
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-3.0000000€ 00 


1.0000000€ 00 


1. 0000000E 00 


.. 





// (standard OS JOB card) 
// EXECALINCON 
EXSUIENK . SYSINADD.* 
ESENCLUDEASYSLIB(üSTVAR) 
pk 

//GO.SYSIN.ADD.~* 


STVAR FOR LUEN TEST 94 


0.0 1.0 DO OO 
0. 0 0.0 SO O 
0.0 080 0.0 10 
Ur. O 5-0 ES =9.0 
00 0.0 S0 0 
20.0 10.0 0.0 ERR 
(blank card) 

F 

2.0 

Bn, 

TO 20 


(blank card) 
/* 
Results shown in Fia. 3-28 indicate that the system is com- 
pletely controllable, the plant eigenvalues are -3, -5, -l, 
and 0, the feedback coefficients are 20, 13, -6 and -3 and 
the controller gain K equals unity. 

Step 2 
The observability index is determined using the subprogram 


OBSERV. The computer cards are as follows: 
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// (standard OS JOB card) 
// EXEC, LINCON 
//LINK.SYSINADDA* 


nn INCLUDE SYSLIB(OBSERV) 


Dex 
GO SYSIN,DD,* 
LUEN TEST 0402 
0.0 1.0 0.0 0.0 
0.0 0.0 deep) 0.0 
0.0 0.0 0.0 1720 
0.0 t0 223-0 =9.0 
i. 0 0.0 0.0 0.0 
0.0 1.0 0.0 0.0 
/* 

DEEp..3 


An observability index r = 3 (results taken from OBSERV 
output, Fig 3-29) permits us to design an observer of order 
equal to or greater than (r-1) = 2. Here a reduced-order 
observer is being designed and eigenvalues of -3.5 and -4.0 
were selected for the observer. Note that, as required, 
there are no common eigenvalues for the plant and the 
observer. 
Step 4 

The data for the subprogram LUEN are: 

system order: 04 

number of measurements: 02 


order of the observer: 02 
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0 1 0 
A = ; plant matrix 
0 0 l 
-15 -23 -9 
T 1 : dÉ A 
b = [0 0 0 A distribution. matrix 
0 0 | 
C = , state measurement 
~ Di 0 0 matrix 
e o 13 z6 -3] , desired feedback 
E coefficients 
(from STVAR) 
observer eigenvalues: -3.5, -4.0 


The following set of cards is then: 
// (standard OS JOB card) 

/ / EXEC ,LINCON 

//LINK.SYSIN,DD,* 


EXINCEUDE „SYSLIB (LUEN) 


/* 
//GO.SYSIN,DD,* 
LUEN TEST 040202 


0.0 SAO =23.0  -9.0 
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43 
From the results shown in Fig 3-30, the complete system can 


be described as: 





1 0 0 x, (t) 0 
2 0 1 0 x, (È) 0 
x(t) = + ` u(t) 
b 0 0 jt x. (t) 0 
3 
-15 -23 -9 x, (t) 1 
K -7.5 -1 x, (t) 
x(t) = . 
5 -14 0 x, (t) 
85.5 29, 25 x, (t = 3 
+ + E) 
0 0 x» (t) -1.5 
Et) = [1.0]r(t) - [20 8.5] xi (t) 
x, (t) 
Mo 0.0] a 
x, (t) 
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RSERVABIL ITY INDEX UCULATICN 
ROBLEM IDEN IFICATICHE CeSERV FOR LUEN TEST 


sesti titre 


THE A MATRIX 


0.0 1.000C000E 00 0.0 0.0 
0.0 0.0 1.0000000E 00 0.0. 
0.0 0.0 0.0 120000000E 00 
0.0 1.5000C00€ 01 -2.3000000E 01 -9.0000000€ 00 
THE C MATRIX 
1.0000000E 00 0.0 0.0 0.0 
0.0 1.C00000CE 00 0.0 0.0 
$999490259$324$495953495522249$2252$ $599 $20992$929922*2$ v 
OBSERVABILITY INDEX 3 | 
Figure 3-29  OBSERV for LUEN Test 
LUENBERGER OBSERVER DESIGN PROGRAM 
PROSLEM IDENTIFICATIUN- LUEN TEST 
$:13):$$ 6$ 9$ 36$93$:3$ $89$ $33 RACIAL AICA it 
ThE A MATRIX 
0.0 1.0000000E 00 0.0 0.0 
0.0 0.0 1-0000000E 00 0.0 . 
0.0 0.0 0. 1.0000000€ 00 
0.0 1.5000C00E 01 -2.3000000€ 01 -9.0900000€ 00 
THE B MATRIX 
0.0 0.0 0.0 1.0000000€ 00 
THE C MATRIX 
1.0000000E 09 0.0 0.0 0.0 ` 
0.0 _1.0000000€ 00 0.0 0.0 


DESIRED FEEDBACK CCEFFICIENTS 
2.0000009€ 01 1.3C00000€ Ol 
OBSERVER EIGENVALUES 


=-6.0000000E 00 -3.0000000E 00 


REAL PART IMAG PART 
-3.5000000E 00 0.0 
-4,0002000E 00 0.0 
COEFFICIENTS OF OBSERVER CHARACTERISTIC POLYNOMIAL 
-IN ASCENOING POWERS OF S 


1.4000000€ 01 


7.5000000€ 00 


SESSASHAASTSS RAI 


THE F MATRIX 


-7.5000000E 00 
-1.4000000E Ol 


1,0C00000E 00 
0.0 


1.0000000€ 00 


THE GI MATRIX 
8.5499969E 01 2.9249939€ 01 
0.0 l 0 


THE G2 MATRIX 


-2.9999971E 00 
-1.5000010E 00 


OUTPUT FEEDBACK COEFFICIENTS 
2.0000000E 01 8.5000038E 00 

CCMPENSATOR FEEDBACK COEFFICIENTS 
1.0000000E 00 0.0 


Figure 3-30 Luenberger Observer Design - Computer Results 
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From these equations a block diagram or a signal flow graph 
could be drawn showing the compensated system. 

Step 5 
The system is simulated by the use of the graphical time 
response subprogram (GTRESP) for a unit step input. After 


some matrix manipulation, the following augmented system is 


obtained: 
0 1 0 0 0 o 
0 0 1 0 0 0} |x, (t) 
E) 0 0 0 jl 0 io 
x(t) -20 -23.5 -23 -9 SE REESEN 
145.5 54.75 0 0 -4.5 1] |x,(t) 
30 27500 0 10,580 x4 (t) 
0 
0 
0 
+ et} 
i. 
23 
uH 55 
y_(t) = [20 10 0 0 0 0 0] x(t) 
x(t) 
m - 9 


l55 











PE 1.0 

ENS step input, r(t) = 1.0 and initial condition 
the initial and final times are 0 and 10, 
respectively, the integration time step is 0.0025, and the 


Porting parameter FREQ is 100. 


The computer deck for GTRESP is then: 
// (standard OS JOB card) ,TIME=2 
// EXEC _LINCONF 
//FORT.SYSIN-DD-* 
SUBROUTINE RFIND(T,R) 
R=1.0 
RETURN 
END 
/* 
7 GENK .SYSINADD.~ * 
ESENCLUDEASYSLIB(GTRESP) 
pex 
//GO.SYSIN.DD.* 


BENEESP FOR LUEN TEST 06 


0.0 I 0.0 0.0 0.0 0.0 
0.0 0.0 1.0 0.0 0.0 0.0 
9.0 0.0 0.0 1.0 0.0 0.0 
=20.0 =23.> 223.0) =o 0 O 0.0 
145.5 54.75 0.0 0.0 -4.5 1-0 
50.0 12275 0.0 0.0 2 0.0 
0.0 0.0 0.0 1.0 23: Zi 
20. 10. 0.0 0.0 0.0 0.0 
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00 1020 0.002 100. 
ES 

The results are shown in Fig 3-31. The user is reminded 
that the observer does supply estimates of the missing 
components of the state vector but at the expense of adding 
its own poles to the over-all system. 

For comparison, a run is also made simulating 
the system that would have been obtained if all states were 
measured, using the feedback coefficients and controller 
gain from STVAR subprogram results. Since the same forcing 
input is used, the control cards remain the same and the data 


cards are changed to read: 


ALL STATES MEASURED 04 


0.0 1720 0.0 0.0 

0.0 0.0 1.0 0.0 

0.0 0.0 0.0 120 

0.0 -15.0 -23.0 = 90 
0.0 0.0 0.0 10 

20.0 20.0 0.0 0.0 

20.0 13.20 -6.0 -3.0 
IO 

(blank card) 

0.0 10.0 0.002 100. 
X 

/* 
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The time response obtained in Fig 3-32 is almost 
identical to the one of Fig 3-31, showing that the observer 
designed does a very good job. 

6. Series Compensator (SERCOM) 

This subprogram is used to design optimal linear, 
time-invariant control systems with incomplete state measure- 
ments. The optimality criterion here is in terms of a 
specified closed-loop transfer function to be achieved. The 
main idea behind the subprogram is to construct a series 
compensator such that the need for feedback from the unmeasured 
state variables is eliminated. The way to accomplish this 
is presented in [8] and [1] and the theory is not repeated 
here. The user should, however, familiarize himself with the 
subject before attempting to solve problems by the use of 
the subprogram SERCOM. 

The following paragraphs outline the computer-aided 
design procedure, the inputs required and the expected output. 
To illustrate the technique an example problem is worked out 
in detail. Notice that the overall procedure differs from 
the one presented in [1]. 

a. Design Procedure 

Before the step-by-step design procedure is out- 
lined, it is necessary to recall the main equations from [8] 
and [1]. First the uncompensated system state equations are 


(as for LUEN) of the form 


159 





sjuswsınsesw 23235 239TdUOI) LIN YNDegp394 AOJ dSHULO SC-C Sandra 


tr 


a ee rn Gr 5 ae DARA AS 0o A 

. AUS e D e e. e a e e e 

e AT . e. o e . e o . . 

e A* e . . e. . . . e. o 

e. A? "e a ə e e. e. e. $ > 

e A . . . e . . . . o 

. ASS o . e e. e . . e o 

. AC . e e e. . e e. e e 

o KS o . o e e . . e e 

Ta --- --- oe ann or PA a m ===- 00 300°8 

. A . e e a i e e. . . . 

e. a e e. e e . . o e * 

. Au . . . o . . o e . 

e pe e e e ^ e e e . > 

e Ae 9 . a e . o ® . 

. E e. e e . e. . . . . 

e. A? e e e e e > o . e 

^ A O o . e e . e . e = 

Wu er uie o NEC CE ee B: cs x 5c Ee 

e. A . e . . . . . . » . 

e Á e . e . a e . . . © 

D A e e o e e . . ® = 

e A = e e . o ə . . E e 

. A e e . . e. . . e e o 

e AUS e. e. . . e. a . © . 

o A ə . e e e. . e . . 

. A a e . e. . e a . e 

ST 11090 cx O A i A eee DD c SEE a Uu. EN co 

e GE o e. e e e. o . . 9 

e x o > e o . . e $ E 

+ A? a . e e. e e e e e 

. À e e D e. D e. a a D e. 

e A e. D D e. e D e. e. . = 

° A . o . o . . o e o . 

E D e. e e. o . a . e z 

A . o . . e . e. . . o 

BA mmm nny mm mmm mn rn mm A A A nn Bn ne nn fn ne ny ne ee ng oe eee ne 0000 0300°2 
A e e e e $ e . . . > z 

; DI Di à E : : : , . 

É : e e. A e. e e . . . e 

d z e e e e. A e e . . - 

e . E S i S : ^ a R : : 

o o . e = È : 5 : A í : 

> o a e. . . a . e. . . È 


KZ 
= mene oT TR ER en on en nn nn on none onen 0°) 
O z66"1l 10-3%3°5 10-36L*8 10-3591 16-3959  iC-34*5*^S8  lu-d4-^* i0-382°€ 10-361*2 10-369? 1 EFSER 
A iNoainJ 


ILUNAS IIUVIUVA 


3SilUdS 3d W43ASAS 


l60 





x(t) = A x(t) +B u(t) 
Y- (t) = CRI 
Me) = C x(t) 


where 


er 

Q 
ct 
Il 


output variable to be controlled (could 
be one of the measured states or a 
linear combination of them) 


C = output variable vector 


vector of measured components of 
state vector 


he 
ct 
Il 


E = state measurement matrix 


An arbitrary dynamic controller 


z (t) D z(t) * e w(t) 


T 


u(t) (t) 


Il 
En 
UN 


is added to the above system. It is to be noted that z(t) 


are defined as 


zT(c) = [u(t) últ) u(t€ ... uU (ey 


w(t) = SUUS ES 
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and 





The complete system then takes the form 


AO = Azt) bult) 
z(t) = D z(t) + e w(t) 
A 


Kix (t) - k) x(t) - ko z(t)] 


w(t) 


The block diagram representation is shown in Figure 3~-33A. 
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Fig 3-33A Serially Compensated System 
with Complete State Measurements 
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It is clear that this closed-loop system does not solve the 
problem since it uses all the state variables. It is possible, 
however, starting from this system, to eliminate the feedback 
from the unmeasured state variables and this is the purpose 

of the subprogram SERCOM. Thus, given the above control 
system, the computer program accomplishes the necessary 


transformations and outputs the new closed-loop system 


x(t) = A x(t) + b ult) 

y(t) = C x(t) 

v(t) = D v(t) + G y (t) + Ker (t) 
u(t) = £"v(t) di SCH 


fer in block diagram form, as in Figure 3-338), 


with 
A = plant matrix 
b = distribution vector 
D = compensator matrix 
G = major loop feedback coefficient matrix 
e = input vector 
K = input gain (a scalar) 
E = compensator output matrix 
gi = minor loop feedback coefficient matrix 
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Fig 3-33B Serially Compensated System with 
Incomplete State Measurements 
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From theory, such a linear compensator can be designed pro- 
vided the order of the controller is at least (r-1), where 
r is the observability index of (A,C) [8]. 

To summarize the above expos& and give a practi- 
cal means of using the method, a step-by-step design proce- 
dure is presented. After the theory of the series compensa- 
tor method has been assimilated, it should be sufficient to 
just follow these few steps and look at the input format table 
for SERCOM to solve any given problem. 

Step l 
The subprogram OBSERV is used to find the observability index 
il Of (A,C). If the system is observable, the minimum order 
for the compensator is then established as (r-l). 

Beep 2 
The D, e" and e" matrices are selected such that their 


dimensions are 
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1 0 0 0 

0 1 ci 0 0 
D = Ie 

0 0 e + 0 iL 





For instance, for a compensator of order one, 


lee 
EC 1 
D = 0 


while for a compensator order equal to two, 


f = [1 0] 


The augmented system is then written as 
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x(t) A 


^ 


b | 0 x(t) 0 USE) 


D' 2266) 


z (t) 0 


e: 


E — 
He 
— 


This form complies with the format necessary to use the 


subprogram STVAR of step 3. 


E, 
Y (s) 
A desired closed-loop transfer function SEES is specified 


for the augmented system. The order of the combined system 
is (n*r-1). For example, suppose that a third order system 
is to be serially compensated. Its observability index, 
found using OBSERV, is r = 2. Then a fourth order polynomial 
must be chosen to characterize the desired closed-loop 
behavior. 

At this point, all states are assumed to be 
available for measurement and the subprogram STVAR is used to 
obtain the controller gain K and the feedback coefficient 


matrix k E and E It is recalled that k^ contains the 


1 
plant feedback coefficients while E contains those for the 
compensator. 

Step 4 
The compensating elements for the augmented system are com- 
puted and the required matrix transformations accomplished by 


the use of the subprogram SERCOM. The final system takes 


the form 
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x(t) 


vit) = 


u(t) 


y(t) = 


^— 


where all elements are 


Step 5 


A x(t) + b ult) 
Di Ge LUE Ker (t) 


T 

fy(t) * g'y(t) 

C x(t) 

given in the output of SERCOM. 


If desired, the compensated system is simulated using GTRESP. 


As for Luenberger Observers, some simple matrix manipulations 


are required to put the equations into the form 


x(t) 


v(t) 


y, (t) 


x(t) 
A d pst) 
ROM oen s 
C} x(t) 
v(t) 
IUE) 
K 
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The graphical time response subprogram with appropriate 
time specifications is then run. 

ER” Input 

The data deck includes all the parameters defined 

for the augmented system. To avoid any mistake, the user 
should refer to the design procedure for comparison. The 
input data cards start as usual with the problem identifica- 
tion, the order of the plant N, the number of measurements 
M and the compensator order (r-1) or greater. The complete 
system matrices are then presented, one row at a time, in 
the following order: A (NXN), b” (1x5), C (MXN), 


D[(r-1) x (r-1)1, e*[1x (r-1)] and £'[1x (r-1)]. On the final 


T 
2 


controller gain K are presented. For a zero steady-state 


cards, the feedback coefficient matrices k a and k and the 
error to a step input, these would be entered exactly as 
they appeared on the subprogram STVAR output. The following 


input format table summarizes the entries required for SERCOM. 


Entry Input Description Format Columns Used 

1 Problem identification 5A4 1-20, 21-22, 
order of the plant (N « 10), 372 23-24, 25-26 
number of measurements = M, 
compensator dimension = (r-l) 
or greater 

2 Plant matrix A (NxN) (one 2510-3 l-10, 11-20, 
row per card for N < 8; Cl E e et 
one row per two cards for 
N > 8) 
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Entry 


teen R T 
i.e., the problem identification, the A, b j C, D, 


Input Description 


Distribution vector b? (Lx N) 
(one card if N < 8; two 
cards if N > 8) 


State measurement matrix 

C (M*N) (one row per card 
for N < 8; one row per two 
Garas fOr N > 8) 


Compensator Matrix 

Du I) x (r-1)} (one row 
per card for (r-l) « 8; one 
row cards per two cards for 
gel) > 8) 


Input matrix el (1 x (r-1)) 
(one card for (r-1) SE 
two cards for (r-l) > 8) 


Compensator output matrix 
£T[1x (r-1)] (one card 
if (r-l) < 8; two cards 


lesen) > 8) 


Feedback coefficients matrix 
[(k.T kjT] (1xN*r-1) (one 
card if^(N«r-1) < 8; two 
cards if 8 < (N+r-1) < 16; 


three cards if (N+r-1) > 16) 


Controller gain K 


Table XV - Input Format Table for SERCOM 


E Output 


Format 


BIET 


8F10.3 


SERBO 


GENS 


8Fl10.3 


EI 


S EPOR 


Columns Used 
l-10, 11-20, 
21-30, etc. 


UNE 
21-305 etc. 


1-10, 11-20, 
21-30, etc. 


1-10; 11:20, 


22230, etc. 


ETO TIE 
21-7305 ete. 


1-10, 11-20, 
21-30, etc. 


First the information given as input is listed, 


c at 


and SS Se matrices and the controller gain K. Next the 


final compensator system matrix D[(r-1) x (r-1)] is printed 


(the user must be careful not to confuse this matrix with the 


Original augmented system matrix D), followed by the minor 


171 





feedback coefficient matrix gi (1xM) and the major loop 
feedback coefficient matrix G [(r-1) xM). 
d. Examples 

Two design examples are worked out. The first 
one is a simple second order system with only one measured 
state variable. The other is the fourth order system that 
was used to demonstrate Luenberger Observers in the previous 
section. 

(1) Example One 

A design of a feedback system is required 


such that the following controllable dynamical equation 


x, (t) 0 dl x, (t) 0 

E x E u(t) 
X, (t) 0 -1 xX, (t) 1 

y(t) = [1 0] x(t) 


has a time response to a step input approximately the same 
as for a second order system with poles at -l + j. 

Step 1.1 
The observability index for the system can easily be found, 
by hand or by the use of the subprogram OBSERV, to be r = 2. 
Thus a first order compensator is sufficient. 


Step 1.2 


The D, Ee and Be matrices are selected such that: 


US 





D = 0 
E 22: 
y 


and the augmented 


system takes the form 


x (t) 0 1 0 x(t) 0 
Bc - NS u(t) 
Z(t) 9 1 l 2 (t) E 
0 0 0 1 
Step 1.3 


The degree of the characteristic polynomial is then three. 
Since the desired response was specified to be similar to a 
second order system with closed-loop poles at -1-4 j and -1-j, 
it seems appropriate to select these roots plus a third real 
root with a large negative value. The subprogram STVAR is 
then used to calculate the required feedback coefficients and 
the gain for roots at -10, -l+j, -l-j. The computer deck 
for STVAR is 

// (standard OS JOB card) 

// EXEC ,LINCON 

//LINK.SYSIN,DD,* 

anA INCLUDE, SYSLIB(STVAR) 

/* 


1/60. SYSIN.DD,* 
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STVAR FOR SERCOMI 03 


0.0 1.0 0.0 
0.0 31.0 0.0 
0.0 0.0 0.0 
0.0 0.0 1.8 
1.0 0.0 0.0 


(blank card) 
E 
10. 
L. dë 
(blank card) 
/* 
The results are shown in Fig. 3-34. 
Step 1.4 
Sufficient information is now available to run the subprogram 


SERCOM. We put together the data: 


order of the plant = 02 
number of measured states = Ol 
compensator order = Ol 
plant matrix De 0 l 

WE 
distribution vector b = [0 1] 
State measurement matrix C = [1 0] 
compensator matrix D = 0 
input matrix e” = 1 
compensator output matrix g^ = L 
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STATE VARIABLE FEECUACK 
PROBLEM IDENTIFICATION - STVAR FOR SERCOM 1 
9.9*492599299995»929*$992554*9$999595922 9299229599» *$f*fÀ4à 


THE A MATRIX 


1 


0.0 1.0C000CCE 00 0.0 
0.0 -1.0090000€ 00 1.0000000€ 00 
C.0 0.0 0.0 
TE B MATRIX 
0.0 0.0 1.0000000€ 00 
TER PLRRLLEREIRLLURTTERBEFTTROT DR HF AN EHER EFT HER 
OPEN-LOOP CALCULAT ICNS \ 
OENOMINATOR COEFFICIENTS - IN ASCENOING POWERS OF S 
6.0 0.0 1.0000000F 00 1.0000000€ 00 
THE ROOTS ARE REAL PART IMAGINARY PART 
-1.CC00000E 00 0.0 
0.0 0.0 
GË 0.0 
THE C MATRIX sosta 
1.C000000£ UO 0.0 | . 0.0 


NUMERATOR CCEFFICIENTS - IN ASCENDING PCWERS OF S 
1.0000000€ 00 l 
BEE RIEF TRDRERSE TER RERTERRE KREDIT DL LIT R RG N 
CLOSED-LIOP CALLULATIONS 
KEY = F suore P 
The NUMERA TOR OF H-EQUIVALERT - IN ASCENOIMG POWERS CF S 


1.0002000E 00 1.0999994E 00 §-4999995E-01 
THE ROOTS AFE REAL PAPT IMAGINAPY PART 
-GoC4U9L35E-01 0.0 


THE FEEDBACK COEFFICIENTS 
1 000990VE 00 5.045935655ÉE-01 5. 4899995F-01 
THE GAIN a 2.000000U€ Ol 
TH€ CLOSED-LCOP CHARACTERTSTIC PDLY'IOMIAL ~ IN ASCENDING POWERS OF S 


2.c0U000JE Ol 2+0C093C000F Ol 1.2000000F ol 1.0000000€ 00 


THE RUCTS ARE REAL PART IMAGINARY PART 


-1.00900007 o 
-1.00000J0E -1.0000000E 00 
1.0000000E 00 


090 
OO» 


-1.C0000uvE 
MAXIMUM NORMALIZED ERRCR 2 6.94E-07 


Figure 3-34 STVAR Results for SERCOM Test One 
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> ko] = D 0.55 0.55], (from STVAR 
output) 


K = 20 (from STVAR output) 
So the control deck and data cards for SERCOM are: 
// (standard OS JOB card) 
WN TA EXEC A LINCON 
WILENK.SYSIN,DD,* 
» INCLUDE, SYSLIB (SERCOM) 
/* 
HAGO .SYSIN,DD,* 
SERCOM TEST ONE 020101 


0.0 TAO 


ib. 0 

1.0 UNO 5 0755 

po. 0 

/* 

From the results reproduced in Fig. 3-35, it is easy to 


determine the final system as 
xt) = fo 1] x(t) + [o] ult) 
0 -1l l 
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BRASCEN TOENTIFICATIONS ^ "'OPRC« resT à 


$92€926294$53-85994252*«t$9$9956559»5t$k22Su**u«ut£sth 
TRE A MATRIX 


0.0 1 
0.0 | 


THE 9 MATRIX 
0.0 1.0C00COCE 00 
THE C MATRIX 
1.0090900F 00 0.0 
THE D MATRIX 
0.0 
THE E MATRIX 
1.0000000£ 00 
THE F MATRIX 
1.0090000E 00 
DESIRSO FEEDBACK CCEFFICIENTS 
1.0009090E 00 5.4995555E-01 5. 4999995E-01. 
THE GAIN = 2.0000000€ 01 
sE AAA KL BE EHRE KU LER DU LS U TA KL RG HR FR A HR 
THE COMPENSATOR SYSTEM MATRIX 
-1.0999S93E 01 
MINOR LOJ? FEEDBACK COEFFICIENTS 
-1.1000000€ C1 
MAJOR LOOP FEEDBACK COEFFICIENTS 
1.0100003E 02 


ARESSESSARAATEAT SE SIL AS ATH SESKSHKICASTHA KAA SE RARE 


eoo 
oo 


Figure 3-35 Serial Compensator Design - Test One 
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45 ' Ihr Ain n ' 


v(t) 


-llv(t) + lOly(t) + (20) (1)r(t) 


Mis) — v(t) - ily(t) 
y(t) = x, (t) 
or, equivalently, 
x(t) = x(t) 
E DE SO 
IÓ O OE) 


Bicep. 1.5 


This last set of equations can be readily used in the subpro- 
gram GTRESP to simulate the system forced by a unit step 


input. From the above equations one gets: 


0 1 0 

A, = [jui -1 1 

101 Den 
Bo 0 20] 
co m [1 0 0] 
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K = 1.0 

x(t) = 0 

to = 0.0 te = 10.0 
dt = 0.002 FREQ = 100. 


We assemble the computer card deck as follows: 
MN Standard OS JOB card) ,TIME=2 
// , EXEC, LINCONF 
//FORT.SYSIN,DD,* 
SUBROUTINE RFIND(T,R) 
AO 
RETURN 
END 
//LINK.SYSIN, DD, * 
xx INCLUDE, SYSLIB(GTRESP) 
~, ENTRY, GTRESP 
/* 
//GO.SYSIN, DD, * 
GTRESP FOR SERCOM1 03 
0.0 1.40 0.0 


FT. O 1.0 
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0.0 0.0 20.0 

1.0 0.0 0.0 

0.0 0.0 0.0 

1.0 

0.0 0.0 0.0 

0.0 10.0 0.002 100. 
Y 

/* 


The time response shown in Fig. 3-36 can be easily compared 
with the actual feedback system where both state variables 
are available (by the use of STVAR and GTRESP) and a decision 
made regarding the suitability of the compensated system. 

Here it is important to note that the method 
increases the order of the system and adds undesired 
poies. For this reason it is always wise to simulate (using 
GTRESP). Another good way to investigate the results is to 
run the subprogram STVAR in open-loop mode for the same set 
of ecuations as for GTRESP. This gives the designer a double 
check on the accuracy of the solution and verifies the con- 
trollability. These ideas are demonstrated in the second 
example. 

(2) Example Two 

The same problem presented for the 

Luenberger Observer example is used here, this time with a 


series compensator. The fourth order uncompensated system 
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SYSTEM RESPONSE 


SYMBOL 


VARIABLE 


OUTPUT 


3.13€-01 4.17€-01 


de C8E-01 


0.0 


6.25€-01 T7.30€- 


5.21€-01 


1.04E-01 


01 B8.34F-01 9,38E- a 
see. Mc 


0.0 


p 4 
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A 


GTRESP for SERCOM Test One 


Figure 3-36 





x(t) = 0 1 0 0 See 0 u(t) 


0 0 d 0 0 
0 0 0 l 0 
Lo -15 -23  -$ 


with measurement equation 


MM (2 0 0 — x(t) 


and controlled output 
YEE) = [20 10 0 OI AXE) 


is to be controlled so the overall time response approaches 
the one that would result from feeding back the states, 1f 
they were all measured, for a fourth order system with 
Bosed-loop poles at -2, -2, -l4 j2. 

Step 2.1 
The observability index is found by the use of the subpro- 


gram OBSERV to be r = 3. Thus the compensator order must be 
at least (r-l) = 2. 


ICEP 2.2 


D; KS and a Matrices are selected as follow: 


Pea j0 2), e = go 13, £° = {1 0 
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and the augmented system becomes 





x(t) = 0 1l 0 0 0 XE) + 0 u(t) 
Et) 0 0 Di 0 0 ES 0 
0 0 0 1 0 0 
0 -15 -23 -9 1 0 
0 0 0 0 0 0 
0 0 0 0 0 1 
v t) = [20 10 0 0 0 0] x(t) 
z(t) 
Step 2.3 


Pole placement is usually dictated by some time response 
specifications. The desired response given here suggests 
that four of the closed-loop poles be located at -2, -2 and 
-l + j2. The two other roots are undesired and a rule of 
thumb is to place them to the left of the desired ones. 
Here -3.5 and -4.0 were selected and the subprograms 


STVAR run with the following control and data cards: 


// (standard OS JOB card) 
// ¿EXEC LINCON 
//LINK.SYSIN,DD.* 

~~ INCLUDE, SYSLIB (STVAR) 


pas 
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//GO.SYSIN-DD-* 


STVAR FOR SERCOM 2 06 


0.0 1.0 00 omo 0.0 ra 
0.0 0.0 ED 0.0 0-0 0.0 
0.0 0.0 VO 1.0 0.0 0.0 
0.0 715: 23:0 O HO 0.0 
0.0 0.0 0.0 0.0 0.0 SO 
0.0 0.0 050 0.0 0.0 0.0 
0.0 0.0 0:0 0.0 0.0 10 
20. 10. 0.0 0.0 0.0 0.0 


(blank card) 


F 


(blank card) 
/* 
The output shown in Fig. 3-37 gives the gain and 
feedback coefficients that would be required if all states 
were measured. 
Step 2.4 
Since some of the states are not measurable, the subprogram 
SERCOM is used to transform the oridinal system into the appro- 
priate series compensated system. The information necessary 


to run the subprogram is: 
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00 30000000*6 


00 30000000°1 


00 3000000 


9299909 
0000 m=mO 


TO 3000000€*? 


0*0 


00 3000000 


ocoooo 
o € e € a se 
000-0900 


OMI 3591 WOH 10] SITNSIY HVALS LE-€ eanbty 


ERROR FETTE UT EHE TEE TRETEN 


0°0 00 30000000*2- 
LMVd ASYNIOYNI Loyd 1V3IY 3” S1NDOY Idi 
10 39009090° 1 10 3000000°? 
S 30 Su34O0d INION3ISV NH - S1NI1D14330) NILVYINNN 
0*0 0*0 10 30000000?*1 TO 30C00000°2 
ceste XiH1VW 2 3HL1 
0*0 0*0 
0*0 i 9*3 
: 0*0 0*9 
0*0 10-235666655*5- 
0*0 00 30666555*^*- 
0*0 00 :0C06000* f- 
AuVd AuUVNIOvVHI ldVd 1v34 3uv SL100M 3H1 
00 30000000*1 
TO 300000065 *1 : 0*0 0*0 0*0 


S 4) SYIMUd $9NION3OSV NI - SI1N3121233300 YOAYNIWONIO 
SNOT LVINITVI dOOT-N3dU . 
"»t*extesttterthtteestutwvektéstx»shotsktttetteétott 
0*0 0*0 0*0 0*0 
Xl diVW Y 2H1 


0 
0 
TEE 10 3000000 
K 00 3000000 


OOOO 


0*0 2 
0°0 S 
= 10 30090095 °1- s 
0*0 * 
0*0 2 
o*1 : 


0 
0 
€ 
0 
0 
0 


02272095 
kk ke LK 


00 3000009 
XIMIVN V 3442 

PELA derrate ed 
€ NOJYIS MUS AVALS - NOI1V2IJI AN20] HITNOYA 
WIVNTIIZI ITUWINVA 31VIS 
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OML 3S9L WODNIS 103 SITNS9Y AVALS  (^3UuOO) LE-E€ SANDTA 


90-3LS°T = 30983 OI7LIVH3SCN WAWIXYH 


0*0 00 31vC6110*2- 
0*0 00 1908 P95*1- 
0°0 VO 32599654 °%C = 
0*0 00 i/651000°%- 
00 34166666°1 10-:tcli8555*5- 
00 3% 166666° I- 10-3E 110665 °5- 

LUVoO AYYUIOYnI advd 1933 34V SLOOY 3H1 

00 3C000000° 1 

10 3000005€ *1 10 35966665 *1 20 3086696€*2 20 32066619** ZO 3165561**6 20 39154661 *.2 


S 30 Su3MCd StNIUNJOSV NI - TVIWOA TIA IIISIYILIYYYHI dI01-035S01) AHL 


T9 388*556€*1 £ NIv9 3H1 
10-39862912*€ 10-3>896826 "8 10-32109110*9 00 3611%120°8 TO 3902%+12€6* 72 TO 36100000*? 


S 1N31 2133302 X2v00334 3Hi 


10-3$/90^5**** 00 36€£69*€€”* 1- 
10- 35190999 °4- 00 39f595EF°1- 
^ 00 36696tt^*1T 10-31569?t€€*1- 
; 00 30086€te£5*1- T0-3T569?tRn*1- 
LUVd AdYNIDVHI ¿dvd TWIN 3uv S100U 3Hi 
10-3*882* 12*€ 00 3011/5981 *€ 10 312165€09*1 TO I6895429E° € TO 301€*5149*€ 10 35100000*? 
S 30 SUIMDA INTUNIDS Y NI - LNIIVAINOI-H JO YOLVEIRNAN 3H1 


says 3" AJA 
SNGEAVINDIWI JOCI-23S0V) 
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order of the plant = 04 


number of measured states = 02 
compensator order = 02 
0 d 0 0 
0 0 1 0 
plant matrix A = 
© 0 0 0 1 
lo -15 -23  -$ 
distribution vector bi = [0 0 0 1} 
state measurement matrix C = it 0 0 0 
0 l 0 0 
compensator matrix D = 0 1l 
0 0 
i : T 
input matrix e = [0 T] 
compensator output matrix g? = [l1 0] 
at T 
ik, k, ] =  [20. 253 8.07 2607 29053 320, 


(from STVAR output) 


K = 14. , from STVAR output. 
The computer card deck is then 
// (standard OS JOB card) 
7 / ~EXEC.ALINCON 
//LINK.SYSIN.DD,* 
AAINCLUDEASYSLIB (SERCOM) 
/* 
//GO.SYSINADD.~* 


SERCOM TEST ONE 040202 
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0.0 0.0 1.0 0.0 
0.0 0.0 0.0 1.0 
0.0 =15. 225. y 
0.0 0.0 0.0 (St 
0.0 0.0 0.0 0.0 
0.0 IO 0.0 0.0 
0.0 1.0 

0.0 0.0 

0.0 1.0 

1.0 0.0 

20. 22170 8.071 2607 239333 7,3214 
14. 

/* 


The computer output (Fig. 3-38) gives the compensated system 


x(t) = 0 1 0 o |] x(t) + fol u(t) 
0 0 1 0 0 
0 0 0 1 0 
lo -15 -23 -9 1 
v(t) = 0 Il v(t) + [88.75 D y (t) 
0 0 0 0 
+ TATE) 
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10-36666€12*€ 


OM], 3S9L - UbTSIA dqoaesueduo) [erades gE-E Sandra 


10-36666286°8 


r 


10-36666690 *9 


00 30000000 1 


*w5x*9»9$9$54$€$064$629$9$9€&5*$6959$55099969$:5009990929$9$225$5$9 


$0-3?21211^»?**€ Y0 -395+79090*2 
TO 3091€1L0*2- TO 3*0F?64n*8 


$1N312141309 »2v90334 a( 01 sUCvwW 
00 3»$5t4165* R- 10 32€98105*6- 
$1N312144302 :2Y80334 aC)! YONIW 


OU 2966564 4° b= TO 3€07916?*1- 
C0 399950351 40-3€C€$595*5c*&8 


Xl j4V: W31SAS YOULUSNIIWID JHL 
**»*x$42vxxt9*ttt b -tstt*tuss*$9stoesno?stts»2s29*52€66909$ 

TO 30000009? 1 e NIVO 2H1 

TO 365502€5*2 TO 360006000*2 
S1v31913323) 33480333 9291520 

0*0 00 30C0N000°1 

KXIYIYN 4 3H1 

00 30097000°1 0*0 
X1Y1VN 3 34 


0*0 0° 


i XiVivn O 2H1 


00 316660108 
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0°0 


1O 3000000€ :2- 
00 30000000*1 
0*0 


00 30003000*1 
0*0 


0*0 


TO 3000000S° I- 
0*0 
9*0 

00 30000000* 1 


0*0 

00 3N0C00000° 1 
Xivivn J) JHA 
0*0 
Xt uiVu A HI 
USO 

c'o 

nea 

0°0 

XiUulv« VY S5HI 


*€*9$**€6699€96*94995625609$4$$9€69*4009$«50962v$29€20$$2 


€ 1S31 HOINAS SEET E WI WOdd 
AVAHDIA4 VIISIO dülgeng idk ŠIAS 





u(t) = [1 EE sgol y(t) 
y(t) = L 0 0 0 x(t) 
0 I 0 0 
Sep 2.5 


Again it is relatively straightforward to rearrange the 


equations in an augmented system form suitable for simulation 


using the subprogram GTRESP. 


given here. 


0 1 
0 0 
0 0 
A = 
- A 
88.75 -20.74 
0 0 
T 
Gees" 10 y 
gain = 1 
es= [20 10 0 
= o 
x (tg) - 0 


For completeness the result is 


0 0 0 0 
1 0 0 0 
0 1 0 0 
E =9 1 0 
0 0 0 IL 
0 0 A d 05 
0 0 0 14] 
0 0 0] 
190 





The time specifications are chosen to be 


and the control and data cards for the graphical time response 
subprogram with a unit step input are 
// (standard OS JOB card), TIME=2 
// . EXEC ,LINCONF 
//FORT.SYSIN DD * 
SUBROUTINE RFIND(T,R) 
R= 1.0 
RETURN 
END 
p 
EEDBINK.SYSIN.DD.* 
^^ INCLUDE ,SYSLIB (GTRESP) 
/* 
EES EN D .* 


SERCOM TEST TWO 06 


0.0 1.0 0.0 0.0 0.0 0.0 
DIO 0.0 1.0 0.0 0.0 0.0 
0.0 0.0 0.0 1.0 0.0 0.0 
ET. 20]. EE SE JL ca 0.0 
88.75 -20.74 0.0 0.0 0.0 he @ 
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0.0 10. 0.002 T00: 


Results in Fig. 3-39 are very similar to those 
obtained for the Luenberger Observer "n The response 
can be compared against the original specifications. If unsat- 
isfactory, the designer can redo the problem using different 
pole locations. As mentioned at the end e the previous 
example, it might be good to find out if any mistake was 
made by verifying the location of the closed-loop poles. This 
is easily accomplished by running the subprogram STVAR for 
open-loop calculations for the above augmented system. The 
data deck consists of the problem identification, the system 
Order, A, b and C matrices and two blank cards. The complete 
computer deck is 
// (standard OS JOB card) 
// . EXEC, LINCON 
Y /DENK.SYSIN,DD,* 
„INCLUDE _SYSLIB (STVAR) 
/* 


//GO.SYSIN , DD,* 


192 





00 350*1 


OML 359L NODYIAS 103 dSANLI 6€-£ SANTA 


1E 


e A e e e > e e . e > e 
o A e e e e o e e e o o 
CR: S : S : : : e : : 
o A e e 2 o e e a e e o 
.- A e e e o o o o o e o 
e A e o e e e e e e e e 
2 A e o o e e. o e o o e 
. A e e . e o e oe o a o 
o a rr pan nn pr ron paran ron po = = en eo O 
A o e o e E o o e e o o 

o e e e e e o o > o KL 
e k o e e e e o e e e e 
o A e . e e = e e e a e 
e A e . e e o e e e o e 
us s : : : : x ; e : : 
e e e e LÀ e e e o e e 
>. 1 m e e e e e e e e o 
EE 
e: ; i = - : 
E : : : S : : : : è 
e A o e e e e e e e e e 
e A e e e e e o e e © e 
o A e e e e e e e e e o 
e A o e > e e e e e è e 
gni enn nang on nnn nnn g nnn nnn nnn gn nnn nny enn nnn nnn nn EES 
e ra e e e e e o : E n a 
e Á o e e e 0 e e e e e 
A >. e a e e e. e o o e 
A e e > e o ^ e e e e 
e À e e e e e e e e o ® 
e A e e e . o e e e e o 
e A e. e e e. e e e e o 
EE ami o_o ipo ee We ze e e mm mm Wm e www we gene EE OCC CS ss fb ees rss ess y mm meo ee "es A em m m o m ce Y 
o A o e > e e e K 

e e e e È e A e e o o e e 
o o o e e e A e ar È E i 
: S : 3 : : : è MA > : 
e. e. a e e. o e. o e e. A e 
e e > e e e . e e e A A e 
Y e o a ° e e e e e e A 
ho #7 7-0 0 1-42 42222 0 - -- = --- - - A pp nn = ba en te eA 


T0-3*5**6 10-36t°8 T0-3%f°} 10-362°9 10-392°S 10-361*% 10-:S1t°€ 10-301*2 10-350*1 


A 


TORWAS 


indind 


ZIAVIUVA 


3SNOdS JU WILSAS 


te 


00 


00 


00 


0°0 


300 *1 


300*8 


300*9 


309°% 


300*2 


0*0 


193 





SERCOM TEST 2 06 


0.0 P0 0.0 0.0 0.0 0.0 

0.0 0.0 SO 0.0 0.0 0.0 

0.0 0.0 0.0 1:0 0.0 0.0 

-54. ALO E. x20 160 0.0 

38.75 ED 0.0 0.0 0.0 1.0 

0.0 0.0 0.0 0.0 EE DNE 
0.0 0.0 0.0 0.0 0.0 14.0 
20. 10. 


(blank card) 
(blank card) 
/* 

Note that the gain K is carried inside the b. 
matrix as required by the equations representing the final 
compensated system. Results presented in Fig. 3-40 show that 
the roots are very close to their originally specified 
locations. 

7. Optimal Control/Kalman Filters (RICATI) 
RICATI is a double-precision subprogram used to solve 


the Riccati differential equations 


B(t) = -P(t)A - A P(t) + P(t) BRB P(t) -Q (1) 
and/or 
Pit) = ab(t) + PA™-- B(t)c’R eg + BOB (2) 
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0°0 


TO 30000056€”1 


TO 3000000+*1 


00 3000000 


suor3e[gnopeo doog-uedo 10] AVALS butsn 
OMI 3591 WODANIS 103] SITNSY9Y JO UOT3EOTJTJ9A OPF-E€ SANSBTA 


i 0°0 00 30000002°?- 
LUVd ADNIOVHI ddva Wis 33y S1098 341 
cO 3303390%°! eO CC OC CU 
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to obtain the gain matrix 


G¿(t) = R BP(t) (5) 


or/and 


1 


G.(t) heer (te) (4) 
Equations (1) and (3) pertain to the solution of the 
state-regulator problem while (2) and (4) occur in the con- 
tinuous Kalman filter algorithm. For convenience a brief 
discussion of each subject is included. First the state- 


regulator problem: given a linear, time-invariant system 


[9] 


x(t) a) o 


— 


y(t) 


Il 
UO) 


x(t) 


where u(t) is not constrained, a control law is to be found 


such that the quadratic cost function 


t 
£ 
B= six (t.)P-x(t,)1 + 5 f(x" (t)Qx(t) «u^ (e) Ru (t) lat 
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is minimized. Such an optimal control exists, provided 
that Pe and Q are positive semidefinite and R is positive 
definite, and is given by 


u(t) = -R71 


B P(t) x(t) 2-G (t) x(t) 
where P(t) is the unique solution of the differential Riccati 
equation 


l 


P(t) = -P(t)A - ATP(t) + P(t)BRU B P(t) - Q 


with the boundary condition P(t,) = Pei t. is a specified 
value. The RICATI subprogram is used to determine the control 


qain matrix 


such that the closed-loop system 


x(t) = Axtt) + Butt) 
u(t) = -G_(t) x(t) 
y(t) = C x(t) 


is optimal with respect to the specified performance measure: 


The computer can solve for either or both the transient and 
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the steady-state control gains.* Notice that the gain matrix 
Gt) output by the computer does not include the negative 
sign of the feedback loop. 

For the second type of problem, a continuous Kalman 
filter is to be obtained and the subprogram RICATI is used 
to find the optimal filter gain matrix for the design. Here 
again the user has a choice of getting either or both the 
transient and the steady-state gains.? The problem to be 
solved is to find an optimal filter for a linear, time- 


invariant system [10] 


x(t) 


Il 
d 


x(t) TB w(t) 


z(t) 


Il 
O 


A le) 


where v(t), the measurement noise, is uncorrelated and has 


covariance matrix Q. The random process forcing input w(t) 


rhe conditims sufficient for steady-state control 
to exist are that the system be completely controllable, i.e., 
the matrix [B AB ... AN-1B] be of rank n where n is the 
order of the plant, that no terminal cost be considered in the 
cost function and that A and B be time-invariant. [9] 


seen E conditions for steady-state filter gains to 
exist are [10]: 
(a) the plant must be completely observable 
(b) the plant must be time-invariant, i.e., A, F 
and C are independent of time ur 
(c) the random processes v(t) and w(t) are stationary, 
i.e., R and Q are constant. 
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is also uncorrelated and has covariance matrix R. The 


expected values of the initial states are 
x9 7 E[x(t)] 
The solution is obtained by choosing the filter gain matrix 


Gei = R P(t) 


such that the plant, measurement and Kalman filter are 


x(t) = A x(t) + B w(t) 
z(t) = C x(t) + v(t) 
x(t) = A R(t) + Gy(t) (z(t) - C Sta) 


These equations are also presented in block diagram form in 
Figure 3-41. 
The purpose of the subprogram RICATI is to solve the 


differential Riccati equation 


Bit) = A B(t) + PaT + BOB" - P(O C TR CPt) 


with initial condition 
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P(t.) = P, = El(x(t,) - x) * (x(t) - T 


to calculate the filter gain matrix Ge (t). 


a. Input 


A cormon input format applies to both state- 
regulator and Kalman filter problems. However the matrix 


definitions differ. 


w(t) 





Fig 3-41 Continuous Kalman Filter Block Diagram 


(1) Basic Input 
The input data deck first card contains the 
problem identification, the order of the plant (N < 10), the 


number of control inputs (M < 10) and the number of measured 
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== 


outputs (L < 10). Since these numbers define the dimensions 
of each subsequent matrix, extra care is suggested. Next 
the plant matrix A (Nx N) is entered one row at a time. Simi- 
larly the control matrix B. (Mx N) and the observable output 
matrix C (LxN) are given. The above forms the basic input 
and needs only be included once. 
(2) Control Option Input 

This portion of the data is used when solving 
state-regulator problems. The letter C is printed in the 
first column of the first card to indicate that option con- 
trol is selected. On this same card, if and only if transient 
gains are desired, the user gives the initial time to” 
the final time te and the number of time points of the control 
gain matrix ( NPOINT ) to be printed. If the steady-state 
solution only is desired, the letter C still appears in 
column one but the rest of the card is left blank. 

Next the control weighting matrix R (M x M) 
is entered, followed by the state weighting matrix Q (Nx N). 
If and only if the transient response of the gains was 
requested, by assigning non-zero values to to” te and n 
points, the terminal boundary condition matrix P (te) (N x N) 
is given last. 

ME iter Option Input 

The first card of this portion of the data 

deck indicates a Kalman filter problem by the letter F 


punched in column one. As for the control option input, 


the time interval and number of points of the filter gain 
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matrix transient response to be output are also entered on 
ERBE first card, if and only if the transient response is 
desired. Next, the measurement noise covariance matrix 

R (L x L) and the random input covariance matrix (Q (M xM) 

are entered, one row at a time. Finally, if and only if the 
transient filter gain solution was requested by assigning 


non-zero values to t t, and NPOINT the initial boundary 


om f 
condition matrix P (tq) (Nx N) is given. 

(4) Problem Termination Card 

The user may ask for several different com- 

puter solutions of the same basic problem by stacking the 
control input cards for transient response and the control 
input cards for steady-state solution, or the filter input 
cards for steady-state solution and the filter input cards 
for transient response.  Terminaticn of a given problem is 
indicated by a blank card. As usual, many problems can be 
executed under the same run by placing the complete data 
decks one on top the other. 


The followina input format table summarizes 


the above. 


Entry Input Description Format Columns Used 
l Problem identification, 5A4, l-20, 
peti order of the plant (N Selo) 12% 21-22, 
number of control inputs 12, 23-24, 
(M < 10), 
number of measurements I2 25-26 
(L < 10). 
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Entry 


Basic 


Basic 


Basic 


5 


Gentrol 
Option 


6 


Gemeen 
Option 


> 

Control 

Option 
8 


ff 
NPOINT 
#0 
Control 
Option 


9 


Filter 
Option 


10 


Filter 
Option 


rr 


Filter 
Optien 


Input Description 


Plant matrix A (NxN) (one row 
per card for N < 8; one row 
per two cards for N > 8) 


Distribution matrix BI (M x N) 
(one row per card for N < 8; 
one row per two cards for N > 8) 


Measurement matrix C (LxN) 
(one row per card for N EO 
one row per two cards for N » 8) 


Letter C, 

initial time t,, 

final time t., 

number of points (NPOINT) 


Control weighting matrix R (Mx M) 
(One row per card if M < 8; one 
row per two cards for M > 8) 


State weighting matrix Q (NxN) 
(One zow per card for N < 8; 
one row per two cards for N > 8) 


Terminal boundary matrix P(t,) 
(N x N) (one row per card for 
N 8; one row per two cards 
for N > 8) 


Letter F, 

initial time to: 

femal time tus 

number of points NPOINT 


Measurement noise covariance 
matrix R (Lx L) (one row per 
card for L < 8; one row per 

two cards for L > 8) 


Random input covariance matrix 
Q (MxM) (one row per card for 
M < 8; one row per two cards 
for M > 8) 
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Format 


8E10.0 


8E10.0 


8E10.0 


Al, 
210237 
P0 
ES 


8E10.0 


82.1050 


8E10.0 


Al, 
P10,3, 
Fl0.3, 
Tes 


8E10.0 


8E10.0 


Columns Used 


1-10, 11-20, 
21-30, etc. 


1-10, 11-20, 
21-30, etc. 


1-10, 11-20, 
21-30, etc. 


1, 
IZA 
21-30, 
3492-33 


1-10, 11-20, 
21-30, etc. 


1-10, 11-20, 
21-30, etc. 


1-10, 11-20, 
217230, etc. 


l, 
11-20, 
212305 
3132233 


1-10, 11-20, 
21-30, etc. 


1-10, 11-20, 
21-30, etc. 








Entry Input Description Format Columns Used 


12 Initial boundary value matrix P (tg) SETOSO l-10, 11-20, 
Iff (Nx N) (one row per card for ^" 21-30, Cre. 
N < 8; one row per two cards 
NPOINT = 
for N > 8) 
#0 
Filter 
Option 
3 (blank card) (indicates problem 8E10.0 (blank) 
termination) 


Table XVI - Input Format Table for RICATI 


D. Output 
The problem identification and the A, B? and C 
matrices are listed for reference. Then the option requested 
is indicated and the R, Q and P matrices are printed. Finally, 
"steady-state solution" or "transient response" is printed, 


followed by the gain matrix G, or Gu: 


f 
C. Examples 
Two problems are worked out to illustrate the use 
of this subprogram. 
(1) Example One 


In the first case we assume the plant 


G x) (t) 


u(t) 


x, (t) 


and wish to determine what must the control gains be to 


minimize the performance measure 
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00 


J = " [a, 4x, ^ (€) + SZ) = Ru? (t)] dt 


where the weighting factors are 411 = 420, = 0 and 


122 
R= 50. The control option is used. The elements necessary 


for the data deck are: 


(Note that C is not used in the calculations but must be 


included since the input table requires it.) 


R = 50. 
i | 
0 


Both the steady-state and transient solution are desired. 


oO 
Il 


For the transient part of the problem, to = 0.0, te = 10.0 
and NPOINT = 020 are selected with the initial condition 


P(t 


P 2 = 0. 


The control and data cards are then 
// (standard OS JOB card) 
// ^ EXEC ALINCON 


/ /TsENK.SYSIN ADD «~* 
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AAINCLUDE-SYSLIB (RICATI) 
/* 
E SECH DD. $ 


RICATI CONTROL TEST 020101 


0.0 120 
0.0 0.0 
0.0 120 
1.0 0.0 
È 

50. 

4.0 0.0 
0.0 0.0 
E 0.0 10.0 020 
30. 

4.0 0.0 
0.0 0.0 
0.0 0.0 
0.0 0.0 


(blank card) 
/* 
The solution in Fig. 3-42 shows the requested steady-state 
and transient response. 
(2) Example Two 
The second problem is to find the optimal 


Kalman filter gain matrix for the following system: 


x(t) - x(t) + w(t) 





OPTIMAL CONTROL/KALMAN FILTER PROGRAM 
PROBLEM IDENTIFICATION - RICATI CONTROL TEST 


009 $ 59 TI tI TITTI tire 


THE A NATRIX 


0.0 1.000000000 00 
0.6 0.0 

THE 8 MATRIX 

9.0 1.000000000 00 


THE C MATRIX 
1.090000000 09 0.0 
SHOSET LESH EN ERS SET EES EEEEES OE SEEERESTAEEREHES 


*.* CONTROL OPTION ssa 


THE R MATRIX 
5.230090000 01 
THE Q MATRIX 
wee 00 E 

STEADY STATE SOLUTION . ` 
GATNS 
2¿.828426140-0l 7.521196010-01 

SESHSSSAES HS HHA SAE SESS SEA ETEK HEISE KE SES ES 


*** CONTROL OPTION sen 


THE R MATRIX 
5.000000000 01 
THE Q MATRIX 


4.000000000 00 0.0 

0.0 0. 
INITIAL CONDITIONS 

0.0 0.0 

‚0.0 0.0 


SESSHTST tia tte 
1.206000006 OL 

0.0 9.0 

1.14C000010 O1 


d 

1.403242550-02 5.4723€4400-03 
TIME =, L.08C000C30 C1 | 

GAINS 

5.62395160C-02 4.438778840-02 
TIME =, 1402C000C40 C1 

GAINS 

1.214259300-01 1,4937398900-01 
TIME =, 9,60C0005C0 CO 

GAINS 

1.93456633D-01 3.028165130-01 
une Be 9.000000630 CO 
2.49254 8240-01 4.775998020-01 


Figure 3-42 
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The random input w(t) is white noise with variance 02=510: 


The observed variable is given by 


z(t) = TI MA vE) 
where v(t) is also white noise with variance R = w. 
From the above it is easy to extract the data necessary to 
solve the problem by the use of RICATI. Writing down the 


elements one gets: 


A i B = [0 1] 
d -4k 

C - [1 0| 

R = 10 

Q = 10 


The initial condition matrix P (tg) is chosen to be 


P (to) = 0 0 
0 0 


The time is specified as being to = 0.0 and te = 0.5 and a 
number of points to be output is NPOINT = 10. The computer 


cards to solve both for the transient and steady-state are 
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// (standard OS JOB card) 
i EXEC LINCON 
/7/LINK.SYSIN,DD,* 

AP TNCLUDE SYSLIB(RICATI) 
re 

//GO.SYSIN,DD,* 

EBCATI FILTER TEST 020101 
0.0 1.0 

0.0 -4.6 


0.0 0.1 


F 0.0 DS 010 
0.0000001 
10. 
B 
0.0000001 
10. 
Es 
Results presented in Fig. 3-43 indicate that the algorithm 
used by the computer to find the steady-state gains is not 
good enough for the problem. The transient response final 
values are used as steady state gains. 
8. Discrete Time Kalman Filter (Kalman) 

This double-precision subprogram is used to calculate 
the discrete Kalman filter gain matrix Gi The theory of 
the discrete Kalman filter can be obtained from many textbooks 
and articles and is not reproduced here. For example see [10] 
End [11]. 
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UPTIMAL CUNTRUL /KAL‘4AN FILTER PRIGPAM 
PROBLEM IDENTIFICATICN - RICATI FILTER TEST 
SOHHSSSS ESSE SEH SeSse ae seHenas Hee sesaeteertsene 


THE A MATRIX 


0.0 1.0000COCCD O00 

0.0 -4.60000000C 00 
THE B MATRIX 

0.0 1.000000000-01 
THE C MATRIX 

1.000000000 00 0.0 


9996969$9609209999099595559569:324969 56949996999 *49090 X 


see FILTER OPTION ses 


THE R MATRIX 
1.000000000-07 
THE Q MATRIX 
1.000000000 01 
INITIAL CONDITIGNS 


0.0 0.0 . 
0.0 0.0 
td III 

TIME =, 0.0 

GAINS 

0.0 0.0 
TIME 2, 4. 999998960-02 

GAINS 

1.717204510 Ol 6.858276180 02 
TIME 2, 9.999997910-02 

GA IPS 

42036759960 Ol 72822250930 02 
TIME =, le499959596S0-01 

GAINS 

3.9+0694580 01 8.080184449 02 
TIME 2, 1.999999530-Ct 

GAINS 

42037878530 Ol 82148048770 0? 
TIME =, 2.495999480-01 

Gains 

4.035261820 01 8.142:48560 02 
TIME =, 2.9999 99370-C1 

GAINS 

42025756700 Ol 89142682460 02 
TIME =, 3.49999927€-01 

GAINS 

4 2035728330 Ol 8. 1435472270 02 
TIME =, 3.9999 99 17C-0O1 

GAINS 

4.035731320 Ol 8.143565030 02 
TIME =, 4. «$99999C60-61 

GA ln S = 

: @ 2035731320 Ol 82143563460 02 

TIME =, 4. 99999896 0-01 

GAINS à 

4,035731300 0l 8.143563600 02 


SFORASS carrer ercer FETTE TE — 


ses FILTER OPTION 898 


THE R MATRIX 
1. 000000000-07 
THE Q MATRIX 
1.000000000 Ol 


UNAB T3 FINO STEADY-STATE GAIN 
PLEASE use TRANS [ENT RESPONSE OPT TON 


Figure 3-43 Filter Option Test for RICATI 
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The following block diagram, definitions and equations 
are nonetheless included to summarize the ideas and clarify 


the notation adopted in this discussion. 





*k|k-1 


Fig 3-44 Discrete Kalman Filter Block Diagram 


From the diagram, one gets the discrete time system 


state equation 


E + IT 
A du tali 


and the measurement equation 
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SM 


Each element can be briefly defined and the matrix 


dimensions noted as: 


X state vector (Nx 1) 


~k 

$ : transition matrix (NxN) 
Wy system random input (Lx 1) 
ro: distribution matrix (N x L) 
Z £ measurement vector (Mx 1) 
Ho: observation matrix (M xN) 
ui: measurement noise (Mx 1) 
Gy : gain matrix (N x M) 


The problem is to minimize 
e; 
1 NE S deque 


with respect to Gi Note that Ji is a scalar." 


4 


J. actually is the trace of the cost function 


l 
J o= Eil Kel) OH 7 Seil 


"> 


where J is a (NxN) matrix. 
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"m 
p 
De 





The solution to the problem is 


= T T ER 
Sk SS (1) 
SE Ee > 
P = è P #3 + 0 (3) 
sn ree 3 
where 
le 7 ^ Xk|k-r ® Sele RR icon] d) 
and 
Ek|k-1 ^7 ? Xk-i|k-i =) 
given the initial conditions 
and 
P = E[S x yor 
-0|-1 = O [SL 


The terms associated with the above equations are 
defined as 
(Nx N) matrix of the covariance of error of 


the estimate at k given observations at 
times vp to and including time k. 


Pk|k * 
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SES : (NXN) matrix of the covariance of error of 
the one-step prediction atk given observations 
at times up to and including time (k-1) 


R : (Mx N) covariance matrix of the measurement 
noise 

Q : (NxN) covariance matrix of the random 
input 


The matrix 


9 - TEMO I 
is computed from the parameters e? (Lx N) and E(W,W, 7] CI ES 
The purpose of the subprogram KALMAN is to solve the 
recurrence relations (1), (2) and (3) for a specified number 
of iterations N and print the filter gain matrix G, as a 
denction of k. 
a. Input 
Since many problems are encountered where the 
designer must compensate for time-varying environment by 
letting the covariance of the observation noise be variable, 
it was decided to permit the user to define the R (M x M) 
matrix with an external subroutine. The subprogram KALMAN 
thus is accessible under Mode Two of operation only. Also 
note that the subprogram is double-precision. 
The first input to be entered is the covariance 
of the observation noise via the double-precision subroutine 


RDEF(R,NP,M) performed by the main program where M is the 


order of the matrix. The parameters NP and M are directly 
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available from the main program and must not be assigned any 


value here, i.e., leave them as NP and M. The subroutine 


is then 


SUBROUTINE RDEF (R,NP,M) 


IMPLICIT REAL*8 (A-H, O-2) 


DIMENSION R(20,20) 


FORTRAN statements defining 


the R matrix 
RETURN 


END 


(see example part c.) 


Next the data deck is punched. The problem 


identification, the order of the 


system N, the dimension of 


the random input vector L and the number of outputs M are 


given on the first card and are followed by the (N xN) matrix, the 


T 


na 


matrix and the initial condition 


accordance with the input format 


Entry Input Description 


1 Problem identification, 
of the system (N < 10), 


PT” (LxN) matrix, the E [WW] (LxL) matrix, the H (MxN) 


matrix P (NxN) in 


-0|-1 
table shown below. 


Format Columns Used 


order SA4, 2I2 1-20, 21-22, 
EE 


dimension of the random input 
vector (L < F0), number of 
measurements (M < 10) 

2 $ (NxN) matrix (one row per 8E10.0 1-10, 11-20, 
card for N < 8; one row per 21-30, etc. 
two cards for N > 8) 

3 2 (Lx N) matrix (one row per 8E10.0 1-10, 11-20, 
Card for N < 8; one row per 21-30, etc. 


two cards for N > 8) 
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Entry Input Description Format Columns Used 


T: 
4 E[WW ] (LxL) matrix (one row 8E10.0 1-10, 11-20 
per card for L < 8; one row 21-30, etc. 
per two cards for L > 8) 


5 H (MxN) matrix (one row per 8E10.0 Jag 11-20, 
card for N « 8; one row per 20123505 etc. 
two cards for N > 8) 


6 Number of time points to be 8E10.0 1-10 
performed (NP) 

7 E L- (N x N) matrix (one row 8E10.0 1-29, 11-20, 
pe tard for N < 8; one row 2l-30, etc. 


per two cards for N > 8) 


Table XVII - Input Format Table for KALMAN 


Dr. Output 
The problem identification, the discrete system 
$ matrix, the transpose of the gamma matrix, the E [WW] 
matrix (listed as the W matrix on the printout), the measure- 
ment matrix H, EE EES 
covariance matrix (at NP = 0) and the initial condition matrix 
are listed for reference. Then the filter gain matrix is 
printed as a function of the time index k, from k = 0 to 
k = NP. 
C. Example 
It is desired to estimate position and velocity 
from noisy position measurements only. The system equations 


are 


x, (k+1) i 45 x4 (k) OL 


= + w(k) 
x, (k+1) 0 1 x, (k) xo 
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z(k) = [l 0] x- (k) + v(k) 


x, (k) 


where the perturbation acceleration w(k) has a root mean- 
square magnitude of 2 meters per second. 


From the above information, one can See that 


D I 5 
0 1 
r7 - [.125  .5] 
e [ww] = mean-square magnitude of the 
ce perturbation acceleration 
= 4, 
The matrix Bol=1 is assumed to be 
P - 10 0 
-0|-1 
0 10 


and the covariance of the observation noise is assumed to be 


l 
5) 


4 ton SNP 


(NP) 


4+ (- thong s PES 10 
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The number of time points to be computed is chosen to be 20 
and the following computer card deck is set up: 
// (standard OS JOB card) 
77 ~ EXEGRLINCONF 
//FORT.SYSIN,DD,* 

SUBROUTINE RDEF(R,NP,M) 

IMPLICIT REAL*8 (A-H,O-2) 

DIMENSION R(20,20) 

DO 1 I=1,M 

DO 1 J=1,M 

IF (NP.LE.10) R(I,J)=4.+(-0.5) **NP 


i PEINE. .GT.10) R(I,J)=4. 


RETURN 
END 
/* 
//LINK.SYSIN,DD,* 


nn INCLUDE SYSLIB(KALMAN) 
ENTRY A SALMAN 
/* 


I GO.SYSUN,. DD, * 


KALMAN TEST 020101 
IL 0:5 
0.0 de) 


02.125 025 


4.0 
20 0.0 
10.0 0.0 


218 





DISCRETE TIME KALMAN FILTER PROGRAM 
PROBLEM IDENTIFICATION = KALMAN TEST 


A ZLBRBRIIZZILRZRIEISE-SEPREIRREPIPRIEPETPTTIITITIITITI| 


The PHI MATRIX 


1.C30000000 DO $5.CCCCCCCCO-01 

0.0 1.0000000CO 00 
THE GAMMA MATRIX 

1.252000090-01 5.00C00C000-01 


THE W MATRIX 
4 .00000000C 00 

THE H MATRIX 
1.00009090C 00 0.0 

THE R MATRIX 
5.,000029000C 00 

INITIAL CCNOITICAS 


1.000090900 01 0.0 


3*0 à232$$€2329$39515$55x:13x x 31 52 x 2: 08 3 0b 309 09 30 ACCU 3 0 


K = D 
GAINS 
6.666656570C-01 0.0 
K = 1 
GAINS 
6.215944510-01 5.5875821*0-01 
K = 2 > 
GAINS 
5.5+504895C-01 5.952565580-01 
K = 3 
GAINS 
6 «24958 2400-01 5.28080571C-01 
K x & 
GAINS 
5.706002310-01 4.261785710-01 
K = 5 
GAINS 
5.47%07715C-01 3.825255070-01 
K 2 6 
GAINS 
5.21220326C-01 3.580885810-01 
Ks 7 
GA INS 
5.l121724200-01 3.522912470-01 
K= 8 
GAINS 
5.06408851C-01 3.506715420-01 
K 2 9 
GAINS 
5.05503350C-01 3.51571€890-01 : i 
K = 6 
eine > SAL ST 475790-01 3.517321640-91 
5.05143227D-01 3.51860€460-01 5 Se È ` 
K= 
CAINS e g 451398560 01 3.517310550-21 
5.05300933C-01 3.520437SS0-91 ° a E 
K = 
Ame 1° e 377300-01 3.517321670-01 
i e e Z = 
5.60530)9330-01 3.519618380-01 y e È 
K = 13 GAINS 
GAINS 5.051374380-01 3.51732915C-01 
5.95255693C-01 3.51€521610-01 8 = 
K = 14 GAINS 
GAINS 5.051375610-01 3.51733025C-01 
5.952039230-01 3.51 777892001 
K= 15 d 
GAINS | 
5.051665970-01 3.51743062C-01 


Figure 3-45 Discrete Time Kalman Bilter Test 
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0.0 15020 
y 
The results of this run are shown in Fig. 3-45. 
Discrete Time Linear State Regulator (STREG) 
This double-precision subprogram is used to compute 
the discrete linear regulator feedback gains F(NS-K). The 
discrete linear regulator problem can be stated as [12]: 


given a time-invariant discrete system represented by 


x(k+1) = Ax(k) + Bw(k) 


where the states and controls are unconstrained, find an 


optimal control u*(x(k), k) that minimizes the performance 


index 
J = > x" (NS) H x(NS) 
N-1 
+3 J (x(k) Q x(k) + w'(k) R w(k)] 
Eco E es 5 Ee 
where 
x(k): state vector (Nx 1) 
A : coefficient matrix (N xN) 
B * distribution matrix (N x M) 


W(k): system input vector (Mx 1) 
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J : performance index (scalar) 


H : real symmetric positive semi-definite matrix 
i (N x N) 
Q : real symmetric positive semi-definite matrix 


(N x N) 


R : real symmetric positive definite matrix (M x M) 
NS : fixed integer greater than 0 (number of stages) 
After solving the problem, one realizes that the 


Optimal feedback gains can be evaluated by solving the 


following two equations only: 


l 


F(NS-K) = -[R+B!P(K-1)B] ~~ [B-P(K-1)A] (1) 


P(K) = [A+BF(NS-K)]*P(K-1) [A + BF (NS - K)] (2) 
T 
+ F (NS - K)RF(NS - K) + Q 
where F(NS-K) is the feedback gain matrix and P(0) = H. 
The STREG subprogram determines the F(NS - K) matrix 
for 0 < NS < 999 as K varies from one to NS. It also gives 
the final value of the real symmetric P(K) matrix, i.e., 


P(NS). From these results the user can design the optimal 


discrete system 
x(k+1) = A x(k) + B w(k) 


w(k) = F(k) x(k) 
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where k = NS-K (a block diagram representation of the system 
is shown in fig. 3-46). Note that P(NS) is presented so one 
can also calculate the minimum cost for the NS-stage process 
given some initial state Xo using the relation [17] 


T 
Eo P(NS) 2 


C4 
dë 
| 
Alt 


O,N .O O 


=a! gigghehbhe 
The first input data card consists of the problem 
identification, the A matrix dimension (N < 10) and the number 
of inputs (M « 10). Then the A (NxN), p! (Mx N), H (Nx N), 
Q (NxN) and R (MxM) matrices are presented one row at a time. 
Finally the number of stages (0 < NS < 999) is given. The 


following input format table further describes the required 


data cards. 


Entry Input Description Format Columns Used 

1 Problem identification, 5A4, 212 1-20, 21-22, 
system order (N < 10), 23-24 
number of inputs (M > 10) 

2 A(N x N) matrix (one row per SEO 1-10, 11-20, 
card for N < 8; one row per 2130, etc. 
two cards for N > 8) 

3 BI (Mx N) matrix (one row per 8E10.0 1-10, 11-20, 
card for N < 8; one row per 21-30, etc. 


two cards for N > 8) 


E H (NxN) matrix (one row per SEl10.0 1-10, 11-20, 
card for N « 8; one row per 21230, tes 
two cards for N » 8) 


5 O (Nx N) matrix (one row per 8E10.0 1-10, 11-20, 


card for N < 8; one row per 21-30, rete: 
two cards for N > 8) 
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x È NS-K 





Figure 3-46 Discrete Linear Regulator Block Diagram 
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Bnery Input Description Format Columns Used 
6 R (MxM) matrix (one row per 8E10.0 OO 
card for M < 8; one row per 230 ,rete, 
two cards for M > 8) 
7 Number of stages for the P3 SS 
processi (0<TNS < 999) 
Dem Output 
The problem identification, the discrete system 
A matrix, the transpose of the distribution matrix and the 
H, Q and R matrices are listed for reference. Then the feed- 
back gain matrix F(NS - K) is printed as a function of the 
backward time index (NS - K) for K = 1 to K = NS. Finally, the 
real symmetric P(NS) matrix is given. 


C. Example 


Given the linear discrete system [12] 


x(k+l) = 0.9974 0.0539 2e e 0-0013 w (kK) 
- „1078 1« 1591 00539 
the feedback gain matrix F(NS-K) is to be determined which 


minimizes the performance measure 


e, 2 2 2 
NEL LL LL 09.25 x. (K) + 0.005 Dä (k)] 
d vn 1 2 
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The data are then: 


e 0.9974 
- .1078 

B? = [0.0013 

ae 

Q = 0.25 
0.0 

R = 0.05 


0.0539 
1. T391 


0.0539] 


0205 


Order of the system, N = 2 


Number of inputs, M = l 


For this problem NS is chosen to be 200 and the 


computer control and data cards are 


// (standard OS JOB card) 
V7 -E Seer GINCON 
//LINK.SYSIN. DD. * 

ente LUDESSYSLIB(STREG) 
p^ 


V7 GO; SYSIN DD. * 
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STREG TEST 0201 
039374 DUIS 
=. 4078 121591 


0.0013 0.0539 


0.0 0.0 

0.0 0.0 

0/25 0.0 

0.0 0.05 
0.05 

200 

/* 


The results presented in fig.. 3-47 show that 


F(NS = K) approaches a constant matrix” F as K > 200. 


a 


IN Muttiple-Input Multipie-Outpue Control System 
Decoupling (MIMO) 


This subprogram is used to determine a feedback con- 


control law 


Wit) =. G r(t) - F x(ty 


"Tf a system is completely controllable and time invariant, 
H = 0, and R and Q are constant matrices then [12] 


F(NS-K) > F (a constant matrix) as NS > œ 


226 





DISCRETE LINEAR STATE REGULATCR PROGRAM 
PROBLEM IDENTIFICATICN = STREG TEST 


TrE A MATRIX 


9.,95740230000-01 
-1.07802000C-01 


TRE 8 MATRIX 
1.30009000C-03 
ThE H MATRIX 


0.0 
9.2 


TRE Q MATRIX 
2.500C00000-01 
0.0 


THE R MATRIX 
5.000006000-02 


5.39 C0C0000-02 


0.0 
0.0 


0.0 
5.0CC0CCCCD-02 


K = 199 
GAINS 
0.2 0.0 
K = 198 | 
GAINS 
-6.101251330-0^ -6.264321€70-02 
K = 197 
GAINS 
-6.532448590=03 =1.473544140-01 
n 195 
-1.564841200-02 -2.612811310-01 
K = 195 
GAINS 
=2.€37855110-02 -4.122C€1C50-01 
Bee" 
-3.121841770-02 -6.082104460-01 
K= 193 
GAINS 
-2.4382354540-02 -8.561395580-01 
Keel"? 
-2.385251940-04 -1.155728250 09 
BS DL 
4 . 951294170702 -1.51735€550 00 
EM 9 
1.297456460-01 -1.920526140 00 
Sege d 
2 .417541560-01 -2.3521€5420 00 
O 
3.826599 1590-01 -2.7925EC4SD 00 
ae 
5.452560770-01 -3.220042300 6023 
le 
7.191423140-01 -3.6113CS5EC 00 
Y 
8.532419290-01 -3.5539£2550 00 
TT d 
1.357520300 00 —4.24151€150 02 


Figure 3-47 


4€**5329€49*935*923*39223*X*933»x-23929x2592552925»x $c xxx 


*9€*94235€ES*92Z2€5XX2€95€99255X2233:352S9*52»:2992x*x»&*5925»x*955£4599*t»* 


iE. e 
PEE EC 


er Er ua 
K = 14 7 
GAINS 
-5.522296540-01 -5.969015290 09 
K= 13 
GAINS » 
-5.,522296540-01 -5:965015950 OO 
K= 12 
GAINS 
-5.52229654D-01 -9.96S015CSC 00 
K = 11 
GAINS 
-5.22296540-01 -5.956S015CSD 00 
K= 10 
GAINS 
-5.522296540-01 -5.569015C90 00 
K = 9 
GAINS 
-5.522296540-01 -9,96S9015CS0 00 
K = 8 
GAINS 
-95.522296540-01 ~5.-S6SOL5CSD 09 
K = 7 
GAINS 
-5.52229654C-01 -5.565015C390 00 
K = 6 
GAINS 
-5.522296540-01 -5.965015C50 09 
K = 
GAINS 3 
=5.52229€54D-C1 -9.969015CS0 09 
K = 4 
GAINS 
-5.522296550-01 -9.96S015CS0 CO 
K = 3 
GATNS 
-5.522296540-01 -=5.565015050 00 
K = 2 
GAINS 
-5,522295540-01 -5.56501£C050D 00 
K 2 1 
GAINS 
-5 „52229654D-0O1l -9,.96S01tCs0 OQ 
K = 
GAINS 
-5.522296540-01 -5.565015C50 00 


AND THE P(200) MATRIX IS 


1.652262530 01 
1.017384130 00 


Discrete Linear State Regulator Test 
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for an Nth order system 


x(t) 


A x(t) + B u(t) 


y (t) 


Il 
de 


x(t) 


such that the control system is decoupled, i.e., the ali input 
r, (t) affects only the ss output y; (t). Notice that the 
subprogram applies only if the number of inputs is exactly 
the same as the number of outputs. The computer calculates 
both the feedback gain matrix F and the command input gain 
matrix G. The user only has to feed in the coupled system 
matrices A, B and C and specify, arbitrarily, the desired 
closed-loop poles of each SEH transfer function Y, (S)/R, (s). 
The theory regarding the algorithm used for decoupling 
is not presented here. For this, the reader is referred to 
[1]. Sufficient information is included, however, to illus- 


trate the concepts and to permit easy use of the subprogram. 
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a. input 

The problem identification, the system order (N 
less than or equal to ten) and the number of inputs and outputs 
(M less than Or equal to ten) are given on the first card 
according to the format shown on the input format table for 
MIMO. Then the A matrix (NxN) is entered, followed by the 
BI matrix (MxN) and the C matrix (MxN), one row at a time. 
Note that B is transposed and the number of inputs must equal 
the number of outputs. 

Next the option card is punched. If the option 
is blank, the phase variable form of the decoupled system is 
obtained and the subprogram returns to begin another problem. 
If options P or F are selected, the control law u(t) necessary 
to achieve a decoupled system with closed-loop poles at loca- 
tions specified by the user is determined. If option =F, 
the next cards give the desired poles of Y1(s)/Ri (s). 

According to the convention established before, 
lf option F is selected the real part of a root is entered 
as being positive if it lies in the left-half plane, negative 
if in the right-half plane and only the positive imaginary 
part of a complex pair is given (see p. wo If option P is 
selected, the coefficients of the characteristic polynomial 
of Y, (s)/R, (s) are entered in ascendina order, the coefficient 
of the highest order term always being unity. 

The subprogram then returns to read option P or 
F and the second decoupled subsystem desired closed-loop poles, 


and so on for the M subsvstems. 
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The design of a decoupled system requires two 
separate runs of the subprogram. First the user must determine 
if it is possible to decouple the system and, if so, obtain the 
order of each decoupled subsystem Y. (s)/R; (s). This is done 
Ey leaving the option card blank. The order of the denominator 
polynomial becomes the order of each decoupled suksystem 
which determines the number of poles or the degree of the 
characteristic polynomial to be selected for closed-loop 
calculations. 

Options P or F are selected for the second run 
and the subprogram computes the control law u(t) which 
decouples the system and places the poles at the selected 
locations. The following input format table summarizes the 


pertinent information. 


Entry Input Description Format Columns Used 
E Problem identification, order 5SA4; 2I12 1-10, 2022, 
of the system (N < 10), number 23-24 


of inputs and number of 
outputs (M > 10) 


2 A matrix (NXN) (one row per 8F10.3 1-10, 11-20, 
card for N < 8; one row per 2-30, etc. 
two cards for N > 8) 

3 Br matrix (MXN) (one row 8F10.3 1-10, 11-20, 
Per card for N < 8; one 21-30, etc. 
row per two cards for 
N > 8) 

4 C matrix (MXN) (one row per 8F10.3 l-10, 11-20, 
card for N < 8; one row 21-30, etc. 


per two cards for N > 8) 
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Entry 


option 


6 


(AE f 
option 
=P) 


7 


IEF 
option 
=F) 


Input Description 


blank = analysis only 

P = closed-loop 
polynomial input 

F = closed-loop poles 
input 


Option 


Polynomial coefficients in 
ascending power of s 

(see complete description 
P. 31) 


Roots of characteristic poly- 
nomial (one root per card) 
(see complete description 

P. 32) 


Format 


Al 


SESS 


8F10.3 


Columns Used 


d 


021712207 
212305 etc. 


l-10, 11-20, 
21-30, etc. 


The above information should become clear from the example 


presented in part c. 


for reference. 


tation of each subsystem is printed. 


Output 


The problem identification, A, E 


and C are listed 


Then the decoupled phase variable represen- 


The denominator poly- 


nomial in ascending powers.of sis given first, followed by 


the numerator polynomial both in unfactored and factored form. 


It should be noted that the subprogram outputs the cancelled 


zeros of Y; (s) /R; (s) as well. 


If closed-loop calculations have been requested, 


by letting option equal P or F, each subsystem closed-loop 


polynomial is printed again both in unfactored and factored 
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form. Finally the feedback gain matrix F and the control 
gain matrix G are presented. 
ln terms of the original system, the resulting 


closed-loop decoupled system is 


x(t) = A x(t) + B u(t) 
HEEN 
u(t) = F x(t) + G x(t) 


or, in a form suitable for graphical time response simulation, 


x(t) 


(A +BF) x(t) + BGr(t) 


y(t) C x(t) 

It must be pointed out that not every system can 
be decoupled. If it cannot, the subprogram is interrupted 
and the message "BSTAR IS SINGULAR - THIS SYSTEM CANNOT BE 
DECOUPLED" is printed. It is also possible that a subsystem 
may be uncontrollable. This is indicated as such on the 
gies listing. 

c. Example 
A two-input two-output system [13] is to be 


decoupled both during transient-period and steady-state. 
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The first subsystem must approach a second order response to 
a step input with a natural frequency of 10, a damping factor 
of 0.4 and no Steady-state error. The second subsystem must 
also approach a second order response to step input but with 
a natural frequency of 4, a damping factor of 0.6 and no 
steady-state error. 

The original system is shown both in block diagram 


and signal flow graph form in Figures 3-45 and 3-49. 





u, (t) y, (t) 








Fig 3-48 Multiple-Input Multiple-Output Control 
System (Block Diagram) 
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d 1/s d l/s D 


0,0 O = C S o 5 
2 
-7 
1 1/s T 10% 
U,o 9 e $ e ® 
g 4 
-4.3 


Fig 3-49 Multiple-Input Multiple-Output 
Control System (Signal Flow Graph) 


The state variable and output equations can be directly 


written as 


x, = -12x, + x, + 2X, 
x, = -7x, + u, 

= = 4X4 = 2X, + X4 
A = -4.3x, $ u, 
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From these, the matrices A, Sc and C are seen 


to be: 
A = En il 2 0 
0 2m 0 0 
4 0 27 1 
0 0 0 -4.3 
B? =" (ao 1 0 0 
0 0 0 1 
cC = 1 0 0 0 
0 0 1 0 


A first run of the subprogram can be made to 
verify if it is possible to decouple the system and find the 
order of each subsystem. The computer cards are 
// (standard OS JOB card) 

V7) Exec ~LINCON 

I EINK.SYSIN DD.* 

na INCLUDE L SYSLIB (MIMO) 
/* 


MIGOSSYSIN . DD,* 


MIMO TEST ONE 040202 

— 1.0 2.0 0.0 
0.0 =7.0 00 0.0 
4.0 0.0 E20 1.0 
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0.0 1.0 0.0 0.0 
0.0 0.0 0.0 1.0 
1.0 0.0 0.0 0.0 
0.0 0.0 1.0 0.0 


(blank card) 
/* 

The result shown in Fig. 3-50 reveals that both 
subsystems are second order. The closed-loop pole locations 
can then easily be selected for each subsystem. For the first 
one, a second order response is desired such that d es 10, 

č = 0.4. Thus, 


SE + 2Zu_s + w 2 = s? + 8s + 100 
n n 


For the second subsystem, the desired response requires that 


Wa = 4 and zZ = 0.6, so 


s^ + 260 S + w 2 - = + 4.8s + 16 
n n 


EE 73.2) 


The computer deck is then modified as follows: 
// (standard OS JOB card) 

// ^ EXECALINCON 

//LINK.SYSIN 4 DD.* 

~~. INCLUDE.~SYSLIB (MIMO) 

Ja 
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MULTI-INPUTE MULTI-OLTPUT PRCGAAM 
PECULEM TOENTIFICATICN = MIMI TEST CNF 


$€"239$959»23099*9€*€9*7$94£t*t£2»€4€95*€5€94*922à5$595235929g9g595*5*t**£ 


TEE A MATRIX 


“L.23000E Cl LeICCCCE 00 2.60030= 00 0.0 
e -7.GOCCOE 00 e OU 0.0 
4.CCCOCE 00 0.0 -2.0N000E 00 1.00000€ 00 
J.0 0.0 Ge -4,300C0E 00 
THE B “MATRIX 
0.0 1.00CCOF 00 C.C 0.0 
0.) 0.0 0.0 1.00000€ 00 


THE C MATRIX 


L.CCCOCE 00 0.0 Let 0.0 
3.0 0.0 1,0099205 09 0.0 


| i: nTi4Ch2412 1 ATC: CEARCPFCESCCALRIZEIRECALIAZEOLE. 
DECIUPLEN PHASE VARIABLE REFRESENTATION 
save SUASYSTEU 1 
DENIMINATIF PELYNCMIAL = Th ASCFANING PCWERS OF S$ 
1.90735£-06 —3.81470E-06 1.CCOOOE 00 ` 
NU4ECATOR PGLYNGYIAL = IN ASCENDING POWERS CF S 
1.003Y0F 00 
suse SUBSYSTEM 2 
DENOMINATOR POLYNCNIAL =- IN ASCENDING PCWERS JF $ 
3.0 0.0 1.000008 00 
NUMERATOR POLYNINIAL - IN ASCENDING POWERS OF S$ 
1.CCOOCE 00 


Figure 3-50 Computer Output for MIMO Test One 
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//G0.SYSIN- DD-* 


MIMO TEST TWO 040202 

me. 0 1.0 2.0 0.0 
0.0 =7.0 0.0 0.0 
4.0 0.0 2290 1.0 
0.0 0.0 0.0 -4.3 
0.0 d 0.0 0.0 
0.0 0.0 0.0 1.0 
1.0 0.0 0.0 0.0 
0.0 0.0 1.0 0.0 
E 

100. 8.0 170 

F 

2.4 DI 

/* 


From the results given in Fig. 3-51, the decoupled 


compensated system can be written in terms of the original 


system. 


x(t) = A + BE x(t) + BG x(t) 
= -12 L 2 0 0 0 
Xi) + re) 
-156 4 12. -2 100 0 
4 0 -2 L 0 0 
36.8 -4  -18.4 -2.8 @ (is 
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MULTI-INPUTe PULTI-=CUTPUT PROGRAM 
PROBLEM IDENTIFICATION = MIMO TEST TWO 
SSOP SS ASH S SHLST AH SHEE EES PEKKA fed SESE 


THE A MATRIX 


-1.20000E Ol 1.COCCOE OO 26999005 09 0.0 
@eVUVIOE JO 0.0 =2.00000€ 00 1.00000€ 00 
0.0 0.0 G.0 -4,30000E OO 
THE B MATRIX 

de) 1.0C0C0E 00 C.0 0.0 

0.0 0.0 C.0 1.00000€ 00 
Thc C MATRIX 

1.00300E€ 00 0.0 Col 0.0 

0.0 0.0 1.090090 09 0.0 


SSSSHES SE FASE HESS SHH ARETE SAAKHESAS AES SHS HkKaAS ESE 
DECOUPLED FHASE VARIABLE REFPRESENTATION 
eee SUISYSTEY | 
DENOMINATOR POLYNOMIAL — IN ASCENDING PNWERS OF S 
1.907355-C6 -3.814 I7)E-Co 1.00000° 00 
HUYERATOR PCLYNOMIAL = IN ASCENDING POWERS OF S 
1.00300£ 00 
99*»* SUASYSTEM 2 
DENOMINATOR PCLYNCMIAL - IN ASCENDING PCWERS OF S 
3.0 0.0 1.000208 90 
NUMERATCR FOLYNOMIABL - TN ASCENDING POWERS OF S 
1.90022CF CO 
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CLNSED=-LOOP PELYNCMIAL = IN ASCENDING PCWERS DF S 
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329 P REAL PAPT — IMAGINARY PART 
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€99$9?$$$$9€9€9922306$99$$$92$$95$$9$9*»99?»9* LERZERSZLPEPSELEATI 
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Figure 3-51 Computer Output for MIMO Test Two 
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yy (t) = x, (t) 
yo (t) = x(t) 


For comparison with the actual results given in 
[13], the decoupled compensated system was simulated 
using the subprogram GTRESP, for a unit step input. Note 
that since GTRESP only allows for single-input single-output 
simulation, two runs must be made. The data for the sub- ` 


program is 


rise =, 1,0 
A = -12 1 2 0 

-156 4 12 -2 

4 0 -2 IL 

36.8 -4 -18.4 -2.8 
Da =o 100 0 0] , for the first channel 
Si = [0 0 0 16] , for the second channel 
C = [l 0 0 oO, | for the first channel 
c = [0 0 1 0] , for the second channel 
EE 


240 





e 
TZERO = 0. TF = 10. 
DT = 0.01 FREQ = 20 


The output y y (t) and yY, (t) were to be plotted for 
the first and second run, respectively. The complete com- 
puter deck for GTRESP follows. 

// (standard OS JOB card) ‚TIME=2 
// n EXEC,LINCONF 
A EORT.SYSIN DD,* 
SUBROUTINE RFIND(T,R) 
R=1.0 
RETURN 
END 
/* 
FEINE .SYSIN,DD,* 


A TNCELUDE,SYSLIB (GTRESP) 


Pant RY ~GTRESP 

/* 

ELIGO. SYSINADD.* 

GTRESP MIMO 04 

iz 1.0 2.0 0.0 


24056: 0 4.0 12.0 ZIO 
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55.8 -9.0 me: 4 


0.0 100.0 0.0 


0.0 020 0.01 


/* 


0.0 Note: 


20% 


these be and € 
Matrices are for 
the first channel 
simulation 


For the second channel, the subprogram is run a second time 


changing the be and c matrices appropriately. Figures 3-52A 


and 3-52B show that the response effectively meets the 


Specifications given initially. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


The eighteen subprograms presented in the third chapter 
constitute the actual linear control subroutine library and 
this thesis, the user's manual that goes with it. The sub- 
programs are easy to access and have proven to work well. 
(They were all tested by solving several textbook problems 
and, hopefully, all the "bugs" have been eliminated.) Very 
little programming is necessary so the user can concentrate 
on control problems rather than worry about computational 
details. The LINCON library is indeed a nice tool for 
analysis and design of linear control systems. 

Furthermore, the library can still be easily improved and 
expanded. Any FORTRAN subroutines can be modified and 
replaced or new subroutines added by following the simple 
instructions given in Appendix A. Note that, as was done in 
reference [1], the subprograms were written to handle 
systems of order less than or equal to ten. This should 
take care of most of the problems encountered. If, however, 
it becomes necessary to solve higher order systems using 
these subprograms, remember that it can easily be accom- 
plished by re-dimensioning the arrays of the appropriate 
subroutines and replacing them in the load module (again 
following the procedure given in Appendix A). 

Finally, the following recommendations should be taken 


into consideration 
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ED by the users, 


- always use the proper job control language 
cards (i.e., the ones described in Chapters 
II and III) to access the subprograms. 

(2) by the future LINCON library "programmer", 

- before making any changes, be certain the 
control cards are exactly those required 
for the job. Double checking with a 
consultant is always a good idea. 

- always keep a copy of the subroutines' 
listings and the card decks. It is not 
possible to obtain any listings or card 
decks from a load module. 

- every modification Sora be documented 
(complete with examples) and the information 
distributed to the users. 

- after changes have been implemented and 
tested, a back-up copy of the new LINCON 
data sets should be created to replace the 
one on magnetic tape (as specified in 


Appendix A). 
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APPENDIX A 


The LINCON Data Sets 


All the linear control subprograms described in Chapter 
III were placed in a load module (pre-processed by the 
linkage editor) on Disk 02 of the Naval Postgraduate School 
OS/MVT IBM/360. Procedures were cataloged so anyone can 
easily access the subprograms under OS Batch. The job con- 
trol language cards to be prepared to use the load module 
linear control library are given in Chapter II. 

The following paragraphs now present the actual content 
of the load module and explain the procedures to 

- modify or add members 

- change the data set's expiration date 

- delete the data sets 

- list the member names and check the disk space 

- compress the data sets. 

Also, since a back-up copy of the data sets was created, the 
procedure to restore the load module linear control library 
is given as well. 

However, before any attempt is made to "play" with the 
load module, it is suggested that the programmer familiarize 
himself with the latest computer procedures and the linear 
control subroutine library. References 14 and 15 should 


also be read. 
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In Content. of The Library, 

The linear control subroutines library contains a 
total of fifty-two subroutines. These are given in Table 
A-l indicating what subroutines are used by the subprograms. 
A detailed description of most of the subroutines is presented 
in [1]. Note that several minor changes had to be made to 
the subroutines in order to implement the library. These 
modifications do not, however, change the purpose or the 
efficiency of the programs. Anyone interested in the pro- 
gramming aspects of the library must utilize both the sub- 
routine listings and reference [1]. 

2. Data Sets Utilities 

It is probable that the content of the LINCON sub- 
routines library will have to be modified at one time or 
another. The following paragraphs outline the procedure 
and give the job control language cards necessary to carry 
out the changes. 

a. Data Set Listing 

The following set of control cards is used to 

list the load module library content and the spacesit occupies: 
// (standard OS JOB card) 
//» EXECAPGM=IEHLIST 
//SYSPRINTA DD ~SYSOUT=A 
7/BBLA DD. UNIT=3330,VOL=SER=DISK02,DISP=SHR 
//SYSIN «DD . * 
~ x LISTVTOC , FORMAT, VOL=3330=DISK02, DSNAME=F0718.LINCON 
x x LISTPDS x VOL=3330=DISK02,DSNAME=F0718.LINCON 
75 
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e) 


BASMAT 
CONOBS 
FRESP 
GTRESP(*) 
KALMAN (*) 
LUEN 
MAIN 
MIMO 
OBSERV 
PRFEXP 
PRTLOC 
RICATI (*) 
ROOTS 
RTLOC 
RTRESP 
SENSIT 
SERCOM 
STRGG 
STVAR 


CALCU 
CHREO 
CHREQA 
DET 
DIVP X 

DMULT X x 
FORM X 
GRAPH X 

HERMIT X X X 

LINEQ X X 
MAXI Br 

MPY X 
MULT x XEX 
NORMP 

PADD 

PEXCG 

PFEXP 

PHNOM X 

PMUL 

POLRT x 
PROOT x x x x 
PVAL X X 

RUNGE - X 

SEMBL X X X X XX 

SIMEQ X X SUN 
SIMUL X X X 
SORT X 

SPLIT x x X X 
STMST x x 

SUBP x 

TRESP x 

VEN TEO x 

YDOT x 

Y8VSX x 


NX X 
re 

NN 

S X 

NX dé 

NO 


A o dë Dé Dé 


x 
x 
x 
x 
x 
x 


Table A-1 Subroutines Cross List 


(*) These subprograms were loaded with all their necessary 
subroutines. Each one of them requires an external subroutine 
(see Chapter III). 


(**) The subprogram MAIN is used to call all the sub- 
programs (operation under Mode Three). It requires 450K core. 
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All the subprograms and subroutines names are 
listed in alphabetical order and the space occupied and 
unoccupied given in terms of number of tracks and number 
of cylinders. 

b. Changing Expiration Date 

The expiration date of the subroutine library 
must be changed approximately every six months. The control 
cards used to perform the task are: 

// (standard OS JOB card) 
// ^. EXEC, PGM=CEXPDATE 
WWéSPRENT DD. SYSOUT=A 
//DD1. DD XUNIT=3330,VOL=SER=DISK02,DSIP=OLD, DSN=F0718.LINCON, 
// x LABEL=EXPDT=yyddd 
f 

where yy=year (e.g. 80) 

and ddd=day (e.g. 365) 
The last expiration date given was 80182, i.e. Ol July 1980. 
The computer centre normally sends a reminder listing the 
data sets that are about to expire. 

C. Adding New Members or Replacing Existing Ones 

The following control cards are required to add 

a new member or replace an existing one: 
// (standard OS JOB card) 

// 4. EXEC 4 FORTCL, PARM. LINK="'NCAL,MAP,LIST' 
//FORT,SYSIN DD — * 


Subroutine to be modified or added 


ES 
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//LINK.SYSLMOD, DD_UNIT=3330 ‚VOL=SER=DISKO2,DISP=SHR, 
// ^DSN-F0718.LINCON (member) 
7S 
where "member" is the name of the subroutine to be modified 
or added. Note that the complete set of cards representing 
the subroutine called "member" must be included. Before 
placing the subroutine in the load module library, the com- 
puter compiles it. If any error is found, the linkage is 
not executed and the new subroutine is not placed in the load 
module. The user must carefully check the computer output 
and make sure the message "member now replaced in data set" 
is printed. If not, he must correct any error and redo the 
procedure. Note that a lack of space can also prevent the 
computer from linking to the load module. If this last 
Situation occurs the user should run the "data set listing" 
(part a) control cards to see how much space is available. 
If sufficient space can be allocated, he must run the 
"compressing data sets" control cards (part e) to release 
any unused space in the data sets and then execute the 
addition or replacement. 
d. Removing Data Sets 

It is sometimes necessary to remove undesired 
members from the library (to create space or erase useless 
programs). The following control cards are used to delete 


one or several members from the subroutine library: 
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// (standard OS JOB card) 
// ^ EXEC» PGM=IEHPROGM 
//SXSPRINT ^ DD^ SYSOUT-A 
//DDl^ DD^ UNIT=3330,VOL=SER-DISK02,DISP=SHR 
Z/SYSIN^ DD^ * 
^ ^ SCRATCH ^ VOL=3330=DISK02,PURGE, DSNAME=F0718.LINCON, MEMBER=memberl 
a Aa SCRATCH» VOL=3330=DISK02,PURGE,DSNAME=F0718.LINCON,MEMBER=member 2 
/* 
where memberl and member2 are the subroutines to be erased 
from the module. Here the programmer must be extremely 
careful while using this utility. Mistakes can be very 
costly (from scratching the wrong subroutine to erasing the 
whole subroutine library). For instance, using 
SCRATCH VOL=3330=DISK02, PURGE, DSNAME=F0718.LINCON 

would erase the entire LINCON subroutine library. Be careful. 

Also note that scratching a member does not make 
the space it occupied immediately available. The "compressing 
data sets" utility must be run to release the space (see 
part e). 

e. Compressing Data Sets 

The following control cards are used to free 
unavailable space in the data set: 
// (standard OS JOB card) 
// ^ EXEC ~PGM=IEBCOPY , REGION=100K 
//SYSPRINT- DD» SYSOUT=A 
//DDl^4 DD x UNIT=3330,VOL=SER-DISK02, DSN=F0718.LINCON,DISP=OLD 


//SYSUT3^ DD «UNIT=SYSDA, SPACE= (CYL, (1,1)) ,DISP=(,DELETE) 
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//SYSUT4 DD UNIT=SYSDA, SPACE=(CYL, (1,1)) ,DISP=(,DELETE) 
V/SYSIN ADD — $ 

^ ^ COPY OUTDD=001,INDD=DD1 

/* 

Note that the use of this utility is somewhat 
dangerous since a power failure or a machine check during 
compression will make the data set unaccessible by any 
program [13]. 

2. Back-up Copy 
A back-up copy of the partitioned data sets was made 
by copying them onto the magnetic tape NPS 705, file 01. 
The control cards that were used to create it are: 
// (standard OS JOB card) 
//ONE ^EXEC^ PGM-IEHMOVE,REGION-80K 
//SXSPRINT ^ DD « SYSOUT=A 
//SYSUT1l1 a DD AUNIT-SYSDA,SPACE-(CYL,(3,1)) 
//DDX DD AUNIT=3330,VOL=SER=DISK02,DISP=SHR 
//TAPE DD xUNIT=3400-3,VOL=SER=NPS705,DISP+(,PASS),DCB=DEN=3 
DEISUWSTNC DD. * in SE 72 
COPY ~a PDS+F0718.LINCON , TO=3400-3=(NPS705,1), A 
FROM=3330=DISK02, TODD=TAPE 
/* uM 16 

Since it is possible that the data setsimay be 
lost one way or another, it is imperative to have such a 
back-up copy. To restore the LINCON data sets, the programmer 


must first re-allocate space by running the following Job 


control cards: 
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// (standard OS JOB card) 
//TWO~ EXEC ~.PGM=IEFBR14 
//DDL ^ DD ~.UNIT=3330,VOL=SER=DISK02,DISP=(NEW, KEEP), 
// ^4DSN-F0718.LINCON,LABEL-EXPDT-yyddd,SPACE- (Cyl, (2,1,10)) 
/* 
where yy=expiration year 
ddd=expiration day 
Finally, to restore the data sets one only has 
to run the program given below. 
// (standard OS JOB card) 
//THREE ^ EXEC « PGM=IEHMOVE, REGION=80K 
//SYSPRINT ^a DD ^ SYSOUT-A 
//SYSUT1 ^ DD ^ UNIT-SSYSDA,SPACE-(CYL, (3,1)) 
//DDX ^ DD «x UNIT=3330,VOL=SER=DISK02,DISP=SHR 
//TAPE ^ DD ^ UNIT=3400-3,VOL=SER-NPS705,DISP=(OLD,PASS) DCB=DEN=3 
ASIS IN a DD a * e" 72 
an COPY» PDS=F0718.LINCON,TO=3330=DISKO0O2, A 
FROM=3400-3=(NPS705,1) ,FROMDD=TAPE 


^ 
yak column 16 
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APPENDIX B 


List of the Sources for the Examples of Chapter III 


The following table lists the references from which 


the examples worked out in Chapter II originated. 


Section 


WEB 


ETIC 


PEED 


Example 


Lac 


Reference 


Eveleigh, V.W., Introduction to 


Control System, p. 568 (#3), 


McGraw-Hill, 1972. 


Shinners, S.M., Modern Control 
System - Theory and Application, 
2ndeea.> ape Se? (#7.26), Addison- 


Wesley, 1978. 


Brogan, W.L., Modern Control Theory, 
p. 35 (#2.10), Quantum Publishers, 


1974. 


Ogata, K., Modern Control Engineering, 
pi. 517, Brentice-Hall, 1970 


Ogata, K., Modern Control Engineering, 
p. 275, Prentice-Hall; 1970 


Kirk, D.E., Optimal Control Theory - 


An Introduction, pp. 34-42, Prentice- 
Hall, 1970. 


Ogata, K., Modern Control Engineering, 
p. 797, Prentice-Hall, 1970. 


Kirk, D.E., Optimal Control Theory - 


An Introduction, p. 28, Prentice- 
Hall, 1970. 


Ogata, K., Modern Control Engineering, 
pp. 7/28729, Prentice-Hall, 1970. 


Eveleigh, V.W., Introduction to 


Control System Design, pp. 353-356, 
McGraw-Hill, 1972. 
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Section Example 


Sud 


6.d (1) 


6.d(2) 


ZOL) 


Restena 


ROC 


Reference 


Eveleignh;, V.W., Introduction to 


Control System Design, pp. 357- 
360, McGraw-Hill, 1972. 


Chen, C.T., Introduction to Linear 
Syscem Theory, P- 296, HOLE, 
Rinehart and Winston, 1970. 


Same as 5.d 


Ki kc DET Optimal Control Theory - 
An Introduction, p. 4l, Prentice- 


Hall, 1970. 


Kwakernaak, H. and Sivan, R., 
Linear Optimal Control Systems, 
pp. 347-351, Wilev-Interscience, 1972. 


Parker,S.R., Digital Control 
Systems (Class Notes), 1978. 


Kirk, D.E., Optimal Control Theory - 
An Introduction, Prentice-Hall, 1979. 


Mowrey, J.T., Compensator Optimization 
in Multiple Input Multiple Output 
Control Systems, pp. 26-27, Master's 
Thesis, Naval Postgraduate School, 
Monterey, 1979. 


Table B-l List of References for the Examples 


Worked in Chapter III 
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E). 


In. 


12. 


13. 


14. 


REFERENCES 


Melsa, J.L. and Jones, S.K., Computer Programs for 
Computational Assistance in the Study of Linear Control 


Theory, 2nd ed., McGraw-Hill, 1973. 


IBM System/360 Operating System: Messages and Codes, 
GC28-6631. 


Shinners, S.M., Modern Control System - Theory and 
Application, 2nd ed., Addison-Wesley, 1978. 
Luenberger, D.G., "Observing the State of a Linear System," 


IEEE Transactions on Military Electronics, Vol. Mil-8, 
pp. 74-80, April 1964. 


Schultz, D.G. and Melsa, J.L., State Functions and Linear 
Control Systems, McGraw-Hill, 1967. 


Eveleigh, V.W., Introduction to Control System Design, 
McGraw-Hill, 1972. 


Luenberger, D.G., "OBservers for Multivariable Systems," 
IEEE Transactions on Automatic Control, Vol. AC-11, 
pp. 190-197, April 1966. 


Ferguson, J.D. and Rekasius, Z.V., "Optimal Linear 
Control Systems with Incomplete State Measurements," 
IEEE Transactions on Automatic Control, Vol. AC-14, 
Nor April 1969. 


Athans, M. and Falb, P.F., Optimal Control, McGraw- 
Her 7, 279656. 


Kies. E., Optimal Estimation: An Introduction to the 
Theory and Applications, 1975. 


Parker, S.R., Digital Control System Course (Class 
Notes), 1978 


SEENEN Optimal Control Theory’ An Introduction, 
Prentice-Hall, 1970. 


Mowrey, J.T., Compensator Optimization in Multiple Input 


Multiple Output Control Systems, Master's Thesis, Naval 
Postgraduate School, Monterey, T979. 


Naval Postgraduate School Technical Note No. 0141-05, 
User Libraries and Source Code Editing Under OS, 
by S.D. Raney, Revised July HE 
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15. 


BS. 


der, 


T8. 


I9. 


Brown, G.D., System/370 Job Control Language, Wiley, 
1977. 


Brogan, W.L., Modern Control Theory, Quantum Publishers, 
2970. 


Ogata, K., Modern Control Engineering, Prentice-Hall, 
1970. Ey Sun 


Chen, C.T., Introduction to Linear System Theory, 
Holt, Rinehart and Winston, 1970. 


Kwakernaak , H. and Sivan, R., Linear Optimal Control 
Systems, Wiley-Interscience, 1972. 
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